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A  NUMERICAL  SOLUTION  OF  THE  HEAT  TRANSFER  EQUATION 


Sh.  Ye.  MIkeladze 
1.  Preliminary  Remarks 


\ 


The  Cauchy  problem  of  the  propagation  of  heat  In  a  one-dimensional 
conductor  has  been  studied  In  an  article  by  Kurant,  Fridrikhs,  and 
Levi  [l].  The  authors  show  that  for  a  difference  equation  approxi¬ 
mating  a  heat  transfer  equation  to  converge,  the  time  step  should 
decrease  In  proportion  to  the  square  of  the  spatial  step.  It  was 
shown  that  the  method  may  be  used  in  solving  multidimensional  problems, 
i.e.,  those  with  two  or  more  spatial  coordinates.  They  do  not  study 
the  law  of  error  propagation  nor  do  they  estimate  the  error. 

Even  earlier,  Richardson's  work  [2]  had  appeared  in  which  the 
question  of  approximating  the  one-dimensional  heat  transfer  problem 
with  boundary  and  starting  conditions  was  developed  in  one  particular 
example  with  the  aid  of  a  difference  equation.  In  this  work  the 
temperatures  in  the  initial  layer  are  calculated  by  a  Fourier  series 
and  then  they  are  found  by  means  of  a  recursion  relation  layer  by 
layer.  The  convergence  of  the  computational  process  used  in  the  work 
is  not  demonstrated  nor  is  the  error  calculated. 

Subsequently  It  proved  that  Richardson's  computational  process 
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diverges  [2].  This  fact  is  explained  in  a  quite  elementary  manner 
on  page  22  of  an  interesting  survey  article  by  P.  P.  Yushkov  [3]. 

Just  because  of  the  incomplete  state  of  the  theory  of  approxi¬ 
mating  parabolic-type  linear  equations  by  difference  equations  and 
because  of  the  multitude  of  interesting  problems  awaiting  solution 
we  considered  it  expedient  to  devote  a  few  articles  to  this  theory 
[4-6]. 

The  aim  of  the  present  work  is  on  the  one  hand  to  set  forth  the 
problem  of  approximating  parabolic-type  equations  by  more  complete 
difference  equations  with  positive  coefficients  guaranteeing  that  the 
computations  will  converge  for  one-dimensional  and  multidimensional 
problems  under  sufficiently  general  conditions,  and  on  the  other  to 
point  out  certain  findings  [7-11]  made  more  recently  in  a  less 
general  formulation  and  resulting  from  findings  made  by  us  earlier 

[5]. 

The  main  attention  of  this  paper  is  devoted  to  deriving  the 
general  recursion  relations  which  allow  us  to  pass  from  layer  to  layer 
and  also  to  evaluate  correspondingly  the  error  of  the  solution. 

We  worked  out  a  method  of  constructing  a  system  of  linear 
algebraic  equations  relating  to  each  other  the  values  of  the  unknown 
functions  in  the  nodes  of  any  two  successive  layers.  In  the  one¬ 
dimensional  case  it  leads  to  simple  new  formulas  of  great  accuracy; 
the  method  is  also  applicable  to  two-  and  three-dimensional  cases  . 
which  up  till  now  have  been  excluded  from  examination.  This  applies 
especially  to  the  three-dimensional  case.  It  is  demonstrated  that 
the  derived  system  of  equations  can  be  solved  by  iteration,  the  con¬ 
vergence  of  the  Iterative  process  being  ensured  for  any  intiial 

values  by  appropriate  selection  of  time  and  space  steps  from  the  In¬ 
terval  or  intervals  of  change  in  them.  In  proving  the  convergence 
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of  the  computational  process.  It  is  assumed  that  the  desired  solution 
exists. 

Formulas  are  derived  for  rhombic  networks.  From  these,  in 
particular,  are  obtained  the  formulas  for  rectangular  and  hexagonal 
networks  and  the  solution  errors  are  estimated. 

Highly  accurate  formulas  are  derived  for  boundary  conditions  of 
the  general  type. 

The  finite-difference  equations  and  estimates  to  be  derived  in 
this  work  may  be  used  equally  in  the  Cauchy  problem  and  in  problems 
involving  boundary  and  initial  conditions. 

In  the  following,  for  brevity,  we  concentrate  only  on  the  problems 
of  heat  propagation  with  boundary  and  initial  conditions. 

Throughout  this  article  we  shall  assume  that  we  are  considering 
bounded,  solid,  uniform,  isotropic  bodies  and  will  not  especially 
stipulate  this. 

The  length  of  the  article  does  not  allow  us  to  analyze  all 
possible  cases  to  which  our  arguments  are  applicable.  Those  who 
have  familiarized  themselves  with  our  findings  will  easily  see  how 
they  may  be  used  in  various  cases  not  examined  in  this  article. 

2.  Error  Analysis 

In  this  section  we  shall  investigate  the  errors  in  the  various 
algebraic  analogs  of  differential  equations  of  parabolic  type  and  de¬ 
duce  sufficient  tests  for  convergence  of  the  solution  of  the  boundary 
problem  for  the  difference  analog  to  the  corresponding  solution  of 
the  parabolic  differential  equation. 

We  shall  assume  the  values  of  the  desired  function  to  be  known 
in  the  nodes  of  the  first  few  layers  and  to  be  £  in  number  and  shall 
examine  the  difference  equations  which  permit  us  to  determine  its 
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values  In  the  succeeding  layers,  layer  by  layer. 

In  so  doing,  we  shall  limit  ourselves  to  the  case  of  a  spatial 
variable  x  and  time  t.  Inasmuch  as  the  theory  which  we  are  expounding 
Is  applicable  to  any  number  of  spatial  variables;  but  for  definite¬ 
ness  our  arguments  will  concern  one  variable.  The  extension  to  the 
general  case  presents  no  difficulty  at  all. 

We  shall  begin  with  certain  definitions  and  nomenclature  permitting 
us  to  condense  further  arguments,  and  above  all  we  shall  agree,  with¬ 
out  stipulating  it  each  time,  that  when  t  >  0  there  exists  a  desired 
solution  of  u(x,t)  satisfying  the  given  differential  equation  with 
boundary  and  initial  conditions,  and  that  all  derivatives  of  it  with 
respect  to  t^  and  x,  up  to  those  orders  which  will  be  used  below,  also 
exist  and  are  continuous  in  the  closed  region  Q. 

As  for  the  solutions  of  multidimensional  problems,  we  shall 
also  always  assume  that  in  those  regions  where  these  solutions  are 
considered  the  assumptions  of  the  preceding  paragraph  hold  true. 

Let  us  take  rectangle  Q  in  the  plane  x,t: 

o<*Si.  o</<7, 

where  T  is  the  time  interval  during  which  the  process  is  being 
studied. 

Let  us  now  examine  the  rectangular  network 

*“*(*" . . t}  . t) 

parallel  to  the  axes  of  the  coordinates,  the  sides  of  the  cells  of 
the  network  being  h  and  i  ( along  axes  x  and  t) . 

Let  us  examine  in  the  network  the  relationship 
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a”,+y'1 

+7  «!*»•+/-»  “*•+/-»+—+  2_,  **•'  “,,y)+ 


+ !?<.»+> 


o,i,...,  —  —  /;  s^i 


> 


(1) 


relating  to  each  other  the  values  of  u(x,t)  at  the  points  (network 
nodes)  lying  on  the  segments 

t~jl,  /-(/  +  /  —  i)/  (2) 

with  its  values  at  the  points  (ih,  ( J+s)  l)  of  the  segment 
t  ■  ( J+s)  i. 

In  the  above  relationship  ( l)  for  all  the  values  of  v,  s,  and 
J  we  shall  have 

(vA.  (x+/)D; 

but  the  coefficients  a  and  3  are  functions  of  the  points  defined  in 
rectangle  Q,  while  R^  g+j  is  the  remainder. 

Rejecting  R^^  g+j  we  arrive  at  the  recursion  equation 

,+y'i  U" •+'-*  +  ( 5) 

*(*,»+>->  t+j-i +•••+ r.  "*»  > 

i»  * .  ^ 

permitting  us  to  find  very  easily  the  successively  approximated 
values  of  the  desired  function  in  all  the  nodes  lying  within  Q, 
starting  from  the  values  on  the  initial  segments  derived  from  Eqs.  (2) 
when  J  -  0. 

Let  us  now  investigate  in  what  cases,  when  l  -»  0,  the  approximate 

values  of  U«  will  tend  to  the  exact  values  u,  in  the  nodes 

i,  s+J  i,  s+J 

of  the  network. 
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Using  Formulas  ( l)  and  (3),  we  can  write 

V" 


•+/.«  »+/-i+ — + y]  *«,>  s*.  * 

_ T  _ 


.  "f  ^  /“  °J  - 


(4) 


where 


—  a, 


Is  the  error,  and 


?•***  &,.»>*  o). 


fr.. 


?•>  •+>"  y  *•<  i+)-i  *f* }  ^  *i»>  t+j—t  “H — + }  \  •»!  f 


(5) 


Let  M  and  6  respectively  designate  the  greatest  values  of 
|Ra  s+j |  and  |aa  g+j|  In  the  rectangle  Q;  let  us  now  limit  ourselves 
to  the  assumptions  that  coefficients  a  are  not  negative,  the  sums 
of  0*  g+j  having  positive  values  for  all 

a,, -»>•••>  **•  i 
not  equal  simultaneously  to  0. 

We  shall  examine  the  absolute  values  of  the  errors  in  the  nodes 
of  the  initial  segments  t  =  0,  I,...,  ( s-l) i  and  the  sides  of  rectangle 
Q:  x  «  0,  x  -  L;  we  shall  designate  the  largest  one  of  them  by  e. 

We  shall  designate  by  *s+j  the  upper  bound  of  the  absolute  values 
of  the  errors  which  occur  on  the  layer  ( s+j)  i  as  result  of  rounding 
off  the  values  of  U^  g+j  computed  with  the  aid  of  Formula  (3)  •  We 
shall  show  that  Formula  ( 4)  will  allow  us  to  investigate  the  complete 
error  g+j  determined  by  the  errors  in  the  initial  values,  the 
error  in  the  Formula  Rj  8+j*  and  the  rounding-off  errors. 

Indeed,  the  fraction 
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*•»*+/-» 


«1*.  •+/-» 


•  is  the  average  for  the  errors 

...  •  ,Ji 

since  all  the  a’s  are  positive  and  do  not  vanish  simultaneously. 

Hence,  when  J  ■  0,  the  absolute  value  of  the  fraction  in  which  we 
are  interested  is  not  greater  than  e  .  Therefore  Formula  (4)  on 
layer  t  «  si  leads  to  the  following  estimate  of  the  error 

I5<>  (I  S  (5  -j-  Af-f- 

Setting  «  1  in  (4)  and  repeating  the  above  reasoning,  we  shall 
verify  that  for  errors  in  the  nodes  of  the  segment  t  =  ( s+l)  l  there 
exists  the  estimate! 

15.-  *+tl  S  *3*  +  (i  +  8)  M-r  ,<v«  ®  +  ^»+i 

and,  in  general,  for  the  values  of  the  errors  on  layer  t  *  ( s+l) l  we 
obtain: 

|C„ .+,!  Si  +  d  +  8 + 8’+  •  •  •  +  *9  M  +  8  +•••+**  * 

We  shall  examine  these  cases  separately 

#<:»  8-1,  5>:. 

When  6  <  1  we  have 

+  .  (6> 

. .  •  • 

where  4  is  the  greatest  absolute  value  of  the  errors  occurring  because 
of  rounding-off  the  values  of  the  solution  in  each  node  of  rectangle 

Q.  ....  ’ 
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The  error  e6^+*,  decreasing  as  increases,  has  almost  no  effect 
on  ^  e+j ,  for  sufficiently  large  values  of  It  vanishes  when  l  -*  0. 
The  error  M(  1  —  6jl  depends  on  the  remainder  term  of  the  calculated 
formula  (3)  and  6.  The  value  of  1  —6  may,  in  particular,  also  be 
infinitesimally  small  in  comparison  with  l.  Therefore  if  M  is  an 
infinitesimal  of  higher  order  than  1—6,  then  M( 1  —  6)”1  will  also 
cease  to  affect  the  error  s+J*  beginning  with  a  certain  value  of 
l.  It  remains  to  examine  the  effect  of  $( 1  —  6)“l  on  s+j,  i.e., 
examine  the  sensitivity  of  the  desired  solution  to  rounding-off 
errors.  We  may,  by  holding  h  and  i  fixed,  make  $( 1  —  6)”1  arbitrarily 
small,  since  we  are  In  no  way  restricted  in  our  choice  of  After 
this,  it  is  easy  to  indicate  to  how  many  decimal  places  one  must 
calculate  using  (3)  >  so  that  the  rcundlng-off  error  is  almost 
imperceptible. 


When  6  *  1  we  find 


+  fr). 


Hence  it  follows  that  in  the  whole  region  Q  in  the  case  under  con¬ 
sideration  there  exists  the  error  estimate 


Thus  proceeding  from  the  values  on  the  initial  straight  lines 


_f  ,1- l) 

and  x’  «  0,  x  »  L,  we  shall  be  able  to  calculate  with  the  aid  of 
Formula  (3)  all  the  s+j* 8  ln  succession,  layer  by  layer,  and  in 
so  doing,  if  the  initial  values  are  approximated  with  an  error  of  order 
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t  with  respect  to  l  and  the  error  in  Formula  (3)  has  an  accuracy  of 
jT+1,  then  we  can  always  make  the  sum  e+TMl”1  less  than  the  previously 
prescribed  magnitude  if  l  and  b  decrease  indefinitely  and  simultaneously. 
It  remains  to  take  into  account  the  value  of  Til”1  generated  by  the 
rounding-off  errors.  The  effect  of  this  error,  as  in  the  preceding 
case,  may  be  made  imperceptible,  if  It  is  calculated  with  superfluous 
decimal  places. 

Finally,  for  6  >  1  we  obtain  the  Inequality 

If,  «l  S  f. + 4±t)  *+.»,  ( 8> 

V  6-1  )  8  —  i 

from  which  it  is  evident  that  if  l  tends  to  zero,  then,  in  order  for 

T 

the  computational  process  to  converge,  the  magnitude  of  £1  must  be 
limited  and,  in  addition,  e  and  (M  +  $)(6  —  l)”1  must  be  arbitrarily 
small  [in  the  best  case  they  decrease  equally  rapidly,  i.e.,  they 
have  the  same  orders  of  smallness  with  respect  to  l(or  h),  when  l  and 
h  decrease  simultaneously  within  limits] . 

Thus  the  presence  of  errors  in  the  initial  values  and  the  effect 
of  rounding-off  errors  affect  the  final  result  most  in  the  third 
case.  In  this  case  no  matter  how  small  the  initial  errors  are,  the 
right  side  of  (8)  may  at  times  become  arbitrarily  large. 

Therefore  we  shall  limit  ourselves  in  all  that  follows  basically 
to  the  cases  where  6  ^l) .  The  question  of  decreasing  the  rounding- 
off  error  has  been  examined  above  and,  in  general,  will  not  be  examined 
further.  In  other  words,  we  shall  assume  that  the  rounding-off  errors 
may  be  made  insignificantly  small,  for  all  practical  purposes,  by 
taking  enough  decimal  places.  Consequently,  we  shall  henceforth 
write  out  the  estimates  of  s+j,  and  shall  at  times  omit  the  round¬ 
ing-off  errors  6(  l  —  6)”1  and  T*l”  *. 
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The  limitation  which  the  choice  of  the  spatial  interval  imposes 
on  the  size  of  the  time  interval  for  obtaining  a  convergent  compu¬ 
tational  process  will  be  stated  later. 

5.  General  Linear  Nonhomogeneous  Equation  of  the 
Second  Order  with  Two  Variables 


Let  us  examine  the  differential  equation 


du 


bp-  +  cu+t  [(«><>)•  • 

ox  • 


(9) 


dt  dx* 

and  seek  the  function  u(x,t),  which  within  the  rectangle  Q  (of  the 
preceding  section)  will  satisfy  Eq.  (9)  and  on  the  border  of  Q,  will 
satisfy  the  initial  condition 


u (x, o)  «=9«a) 


(10) 


and  the  boundary  conditions 


“(M=  ft  (I).  «•  (L,  0 =/,(/)• 


(11) 


We  shall  assume,  for  generality,  that  the  coefficients  a,b,c,  and  g 
are  continuous  functions  of  the  point  (x,t)  in  Q;  we  shall  assume  that 
the  functions  <p(x),  fi(t),  f 2( t)  are  also  continuous  in  the 
corresponding  intervals 

In  addition  we  shall  assume  that 

- :  f(o)— A(o)aad  ?(0  =  /»(°)* 

Equation  (9)  may  be  replaced  [5]  by  the  finite-difference 
equation 

U„  »+«<+»’»  (12) 
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Estimates  suitable  for  investigating  convergence  can  be  obtained 

almost  at  the  outset  by  determining  the  coefficients  a  of  Formula 

# 

(12)  in  conformity  with  (4)  and  then  by  finding  M,  3^  s+j  and  6. 
For  the  coefficients  of  a  we  obtain  the  formulas 


5  In. 


'(>» 


h* 


-j-  l  C,%  |i 


la «.»  *<»» 

*«■' — £H,+T-^7> 

T  «4>  )’ 


Here  the  previous  nomenclature  is  kept;  thus,  for  example  ^  is  the 
value  of  a  in  the  node  ( ih,  kl)  . 

Then  we  obtain  for  ( 12) 

M-±  (6lMt+h'Miem..  +  2h'M,  |*| _ ),  •  ( 15) 

where  M 2,  M3,  and  M4  designate,  respectively,  the  maximum  absolute 
values  of  the  derivatives 

0*u  (.r,  l)  <Pu  (x,  1)  d{  u  (x,  /) 

dfl  ’  7x*  ’  dS~7 

in  Q.  Further,  by  using  Formula  ( 5)  we  obtain 

P;**+/==  1  4“  I***  •+#* 

6=1  +/  \c\m~t 
T 

and  estimate  (8)  in  this  case  assumes  the  form 

Ifr-wl  <«"*'-  +  (14) 

where  M  is  given  by  Formula  ( 12) ,  and  e  designates  the  upper  limit 
of  the  absolute  values  of  u(x,  t)  in  the  nodes  of  the  network  lying 
on  the  sides  of  the  rectangle  Q: 
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t  ■  0,  X  -  0,  X  m  L 


To  ensure  the  validitiy  of  the  estimate  obtained  ( l4) ,  i  and  h 
should  be  chosen  so  that 

3/«.» 


i  — 


A* 


.  t  +  ±A'„>0, 

2  «(,* 

,  A  it.t 
I  — - - - J>o, 

*  «<>» 

In  all  the  nodes  lying  within  the  rectangle  Q.  The  last  two  In¬ 
equalities  occur  for  all  values  of  a  and  b  and  small  values  of  h.  If 
we  are  Interested  In  a  value  of  i  for  some  specific  value  of  h,  then 
for  convergence  of  the  computational  process  we  must  take 

A* 


provided  that 


2a,,  *—  A* £*,*><> 


in  all,  the  Internal  nodes  of  Q. 

Now  let  the  coefficients  of  Eq.  ( 9)  ,  a,  b,  and  c,  be  constant 
while  c^  0.  Then 

£>■»*+/“ i* 

and,  consequently,  when  c  <  0  we  are  dealing  with  the  first  case  of 
section  but  when  c  =  0  we  have  the  second  case.  Therefore  if  the 
calculations  are  performed  with  the  aid  of  the  formula 

— |^  —  2 -j- /f)|  l/j,,  +  ^  X -j- + 

+  (* ~ ^*)  +  /^»1* 

then  In  order  to  obtain  an  estimate  suitable  for  any  c  £  0  we  must 
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use  (7),  although  for  c  <  0,  In  general.  Estimate  (6)  is  more  ad¬ 
vantageous. 

Accordingly  we  obtain 

i&..+y|  +  -I«/Af,4-A*4M4  +  aA*|*|A#,+mr,) 

provided  that 


Thus 


A* 

a  a  -J-  A *  |c| 


s*  +  +  •"•+ )*■+ ■ -T-f 

Finally  let  us  examine  the  differential  equation  of  heat  propa' 
gatlon  of  a  thin  thermally  Insulated  rod  of  length  L  (with  coeffi¬ 
cient  of  heat  conductivity  a2) : 


dt  ** 


(15) 


under  Initial  Condition  ( 10)  and  Boundary  Conditions  ( 11)  . 

The  finite-difference  equation  approximating  this  equation  has 
the  form 


*  +  (U^i+ 

\ 


(16) 


The  computational  process  converges  when  l  and  h  satisfy  the  equa¬ 
tion 

.  (17) 

A*  —  a 


The  estimate  of  the  error  in  this  case  has  the  form: 


«»AftA» 

i 


-13- 


4.  Numerical  Solution  of  Heat-Conduction  Equation 

In  this  section  we  shall  examine  an  Improved  difference  method 
of  solving  the  general  differential  equation  of  heat  conduction 

a1  -by  (og  vr;/.;  oS<S7),  (  18) 

dt  dx* 

where  b  is  a  positive  constant  and  u(x,t)  denotes  the  temperature. 

We  arrive  at  Eq.  ( 18)  when  studying  the  distribution  of  heat 

In  a  bounded  rod  (0,  L)  ,  If  we  take  into  consideration  the  transfer 
of  heat  to  the  external  space.  The  case  b  «  0  (i.e.,  the  case  of 
the  thermally  Insulated  column)  leads  to  Differential  Equation  ( 15) . 
Difference  equations  of  high  accuracy  approximating  ( 15)  were  studied 
in  another  work  [5]j  later  a  whole  series  of  works  waa  devoted  to 
them  among  which  we  may  note  some  listed  in  the  references  [7-13]* 
Although  Eq.  ( 18)  is  brought  to  the  form  ( 15)  by  substituting 
u  =  e'^v  into  it.  It  would  nevertheless  be  incorrect  to  neglect  the 
problems  of  direct  approximation  of  the  general  equations  by  using 
difference  equations  of  high  accuracy. 

The  basic  question  which  arises  here  consists  in  finding  the 
conditions  under  which  the  computational  process  converges,  since 
the  method  of  approximating  the  differential  equation  by  a  difference 
equation  remains  the  same  as  before  [5]. 

Convergence  occurs  if,  in  addition  to  the  limitations  of  Section 
2,  the  product  lb  (l  as  before  designates  the  spacing  with  respect  to 
t)  is  changed  in  a  special  segment  ensuring  a  bounded  change  in 
X.  *  h2/l a2,  where  h  is  the  spacing  with  respect  to  x.  An  exact 
enumeration  of  the  conditions  under  which  the  process  converges  is 
given  later.  ' 

In  order  to  derive  an  equation  of  great  accuracy  let  us  expand 


In  accordance 


the  differences  u1+1^  ^  ^  and  _1#  k  -u^  k 

with  Taylor's  formula,  set  up  the  expression  u1+1^  ^  —  2Uj.  jt  +  ui  _ 1,1c* 

and  replace  in  it  the  derivatives  with  respect  to-x  according  to  the 
formulas 


the  lower  one  of  which  is  obtained  by  differentiating  the  tipper  one  and 
by  expressing  the  derivatives  contained  in  the  derived  equation  in 
terms  of  the  right  side  of  the  upper  one  and  the  derivatives  of 
lower  order. 

We  find 

,  (  ,  h'b  .  h'V\ 

”•+>■*  +  *-i.»-(^2  +  -jr + 7n* /*'  *" 

/>»  A*  h*  ’  fu&f) 

*"\«*  *"  6a*  J  dt  12 a*  01*  *  360  <*{• 

(<  -  1) A <5  <(«  +  !)*• 


Let  us  now  expand  the  differences  u^^  ^*4  ”  UI  k  ui  k+2  — 

—  u.  .  according  to  the  same  Taylor  formula,  multiply  them  respect- 
±j  K 

ivcly  by  the  arbitrary,  as  yet  undetermined  factors  a  and  P,  and  add. 
We  obtain 

* »+t -j-  j3  v„  »+1  —  (*  +  ?)  M<>  *■=(* + 2P)  l  **4- 

of 

+(.T+1?J'“5?-+—  — 3?— 1 - 5p— 


Subtracting  the  two  latter  equations  term  by  term  we  obtain  a 
relationship  permitting  us  to.  write 


where 


.  *  «<*  +  P  Ui,  —  («<+ ,» t  +  *<-»•  0  + 

+( 2  $ + + t^t)  !>•-“  **• “ 


(19) 
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(20) 


ft»J  ■ 


A* 


<i -aft), 


P-0  -ft) 


/a*  6/V 

A*  A* 


ia/V  ala* 
the  remainder  term  assuming  the  form 


«/*  rf’ii  (*,  ij)  (  4  p/1  (*,  C)  A*  d*  *  (g,  i) 


<hj* 


In  addition,  we  find 


&/-?• 


_A*A _ Aft» 

«*  iaa* 


360  dg* 


(21) 


In  order  for  the  computational  process  to  converge,  we  require 
that  the  quantities  a,  p,  P|  ^  and  the  coefficients  of  satisfy 
the  inequalities 


«l r°.  ?>o, 

«+»*.+■* 1*.T 


1 2a 


,« ' 


In  this  case,  PJ  ^  ^  P,  i.e.,  6  1  and  consequently  Estimates 

(6)  and  (7)  will  occur,  respectively. 


We  thus  obtain  inequalities  of  the  form 

6  — (1— A)X<o, 
n-(i-ji)iSot 
*  +  »A-0  -  A  — A»)X<o, 
24  _ (,g  -  i2A)X  +  (1  —  3a  4-  A«) X*  S o. 


(22) 


where 


A -ft. 


(2?) 


Since  we  are  interested  in  the  case  b  ^  0,  we  r'-iould  attribute 

only  non-negative  values  to  A.  On  the  other  hand,  X  should  be  _ . 

bounded.  Therefore  A  should  assume  values  from  the  Interval  bounded 
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by  0  and  the  lowest  positive  root  from  the  number  of  roots  of  the 

following  equations! 

i— aA-o,  i— A-o.  i  —  A-  A'-o,  ,_}A  +  A*-o. 

By  calculating  we  find  that  this  root  is  equaljto  0.3819660. . . , 
and  consequently  l  must  be  such  that 

0^03819660...  (24) 

and  it  is  obvious  that  if  A  and  X  are  chosen  so  that  the  second  of 
the  inequalities  from  the  top  in  (22)  is  satisfied,  the  first  in¬ 
equality  in  (22)  will  also  be  satisfied. 

Thus  the  first  inequality  can  be  omitted  and  the  possibility 
of  the  remaining  three  inequalities  for  A1 3  defined  on  segment 
(24)  may  be  investigated.  Thus  we  should  investigate  .the  functions 
defined  by  the  inequalities 

F,(A,X)=ia-(t-aA)Xgo, 

Ft  (A,  X)=6  +  12A  —  (1  -  A-  A*)  x  <  0, 

F,(A,X)=a4  —  (»8  —  iaA)X-f  (1  -  3A  + A*)*1^0- 


Investigation  shows  that  in  the  region  D  bounded  by  the  segments 


A  —0,3819660. 


and  the  curves 


12 


I  —  2A 

%  9  — 6A  +  ^  57-56A  +  »A* 

A  3 


the  functions  Pi,  F2,  and  Pa  have  negative  signs  and  that 

- - ■  ta 

1 

S 'CP  “S-c 


—  6A  +  1^  57  —  }AA-f- xaA* 
1-2A  1  —  3A  + A* - 


(25) 


-IT- 


for  any  A  of  (24)  and  any  point  (A,  X)  lying  in  region  D. 

When  A  ■  0,  i.e.,  when  b  «  0,  we  hence  obtain,  in  particular, 
that  for  the  computational  process  to  converge  for  f 15)  realized  with 
the  help  of  Formula  (19)  the  values  of  X  -  h2/la 2  should  be  taken 
from  the  interval 

')  -r  Vsi  • 

This  inequality  was  obtained  by  another  method  by  P.  P.  Yushkov  [12]. 

The  estimate  of  the  error  in  approximating  ( 18)  by  the  difference 
equation  [obtained  from  another  article  [19]  after  neglecting  the 
remainder  term]  has  the  form 

!c<,  .+;IS  ~i$r\  ^  '  ~r  ~  ( 26) 

where  M3  and  Me  designate,  respectively,  the  greatest  absolute  values 
of  the  partial  derivatives 

dt *  rf** 

in  the  rectangle  Q,  and  e  is  the  greatest  value  of  the  absolute  error 
in  the  initial  values  of  u(x,  t)  in  the  nodes  of  the  segments 

t  =  0,  t  =  Z,  x  *  0,  x  ■  L. 

Let  us  note,  moreover,  that  the  values  of  u(x,t)  in  the  nodes  of  the 
segments 

t  =  0,  x  *  0,  x  =  L 

are  known  to  us  from  the  boundary  conditions.  The  values  of  u(x,  t) 
in  the  nodes  of  the  segment  t  =  l  may  be  calculated,  for  example, 
by  using  Taylor's  formula  [5 >•  13] •  ; 

Hence  the  following  conclusion: 

Theorem.  If  A  »  lb  and  X  *  ha/l&2,  during  a  change  of  h  and  X,- 

satisfy  Inequalities  (24)  and  (25)*  then  Eq.  ( 18)  with  Initial 

V  >  -:-E  ~  7'  ’  sro?  t’fpF  " 


-18- 


Condition  ( 10)  and  Boundary  Conditions  (ll)  may  be  solved  numerically 
with  the  aid  of  Formula  ( 19)  without  the  remainder  term  k*  The 
error  |(x,  t)  in  each  internal  node  of  rectangle  Q  at  the  moment 
t  ■  T  will  satisfy  Inequality  ( 26)  * 

In  ( 19)  and  ( 20)  let  us  now  assume  0  -  0  and  transform  the  formulas 
’  obtained  by  using  relationships  ( 23)  Then  the  finite-difference 
equation  approximating  differential  Eq.  (8)  is  transformed  into 

( c/t+i, ,  +  Ut-u  *)  +  u<>  *  (A  <0;  (27) 

the  error  in  approximating  this  formula  is 

(»  —  A)*_ _ (l  -A)’  h*  <Pu  (x,  i>)  (I  -  A)»A*  d* « (g,/) 

6  *’*  1296  a*  d?*  2160  dp  * 

so  that 

L  1296a*  2160  j 

By  computation  we  verify  that  for  (27) 

A*  —  A  +  1  >0.  • 

But,  since  0+j  has  a  maximum  at  the  point  A  *  0  and  equals  unity, 
then  6*1.  Thus,  if  we  solve  the  general  heat-conduction  equation 
with  the  aid  of  ( 27)  ,  we  arrive  at  Estimate  (7)*  where  e  is  the 
largest  absolute  error  in  the  initial  data  given  on  the  segments 
t  *  0,  x  —  0,  and  x  *  L. 

The  following  rule  for  the  numerical  solution  of  Eq.  ( 18)  with 

Conditions  ( 1C)  and(ll)  stems  from  the  preceding  arguments! 

— -  •  Rule.  To  solve  Eq.  ( 18)  numerically,  the  values  of  U(x,  t)  are 

calculated  successively,  layer  by  layer,  with  the  aid  of  recursion 

.formula  *  • 

-  l/«,  »+i“— "f*  — Mi+n*  J  3*  •••)•  _  ____  ? 
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where 


6 

i  — A 


h* 

1 * 


(o^SA<  i). 


and  also  the  values  of  U.  .  In  the  nodes  of  the  first  layer  are  found 

J-l  * 

from  the  values  of  c,  in  the  nodes  of  the  zero  (initial)  layer, 
known  from  ( 10) .  For  the  error  which  occurs  from  using  this  recursion 
formula  at  the  moment  t  *  T  the  following  estimate  holds  true 


15  (*,  t)\  s f  +{JL-J^L*  1  h*. 

I  216  a •  560  J 


The  preceding  rule  can  be  formulated  as  follows  when  A  -  0: 

For  a  numerical  solution  of  Eq.  ( 15)  with  Conditions  ( 10)  and 
(ll)  the  values  of  u(x,  t)  are  calculated  in  the  nodes  of  the  network, 
layer  by  layer,  with  the  aid  of  the  formula 


U(x,  t  -h  Q-  + 


(28) 


and  the  inequality 

V  IJ5  1  / 

ensuing  stemming  from  ( 7)  allows  us  to  estimate  the  error  in  the 
numerical  solution  of  Eq.  ( 15)  with  the  aid  of  (28)  in  any  node  Q 
at  any  moment  t  ^  T. 

There  is  one  more  conclusion  from  formula  ( 28) .  We  arrive  at 
this  formula  by  setting  X  =  6  in  ( 16) ,  inequality  ( 17)  being  satisfied 
at  this  value  of  X,  and,  consequently,  the  computational  process 
converges.  But  now  we  can  no  longer  rightfully  assert  that  at  the 
moment  Jb  the  error  resulting  from  rejecting  the  remainder  term  is 
infinitesimally  small  and  has  at  least  the  fourth  order  of  smallness 
with  respect  to  h,  as  follows  from  (29) •  This  is  a  disadvantage  of 
this  conclusion. 

"•  ”  '  rc'rV:\i  . . ~  ~  Ttcp  H£i7  ” 
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When  we  finish  the  present  section,  we  shall  try  to  reduce  the 
description  of  the  derivation  of  the  high-accuracy  difference  formulas 
for  ( 15)  to  the  widest  limits  of  generality. 

Using  the  expansions  u(x,  t)  according  to  Taylor’s  formula 
t  both  with  respect  to  both  x  and  with  respect  to  t,  it  is  possible, 
by  repeating  the  author's  arguments  [5,  p.  83],  to  derive  the  formula 

AxUi,  A,Ui, -j-  AjJit  »+«  —  + 

+  (2  —  At  —  At~r ...  —  Am)  (/(,■  —  Rt, i, 

where  the  coefficients  Ai,  Aa, . . . ,  Ar  satisfy  the  equations 

-f  ..."f 

At  -4"  2*iJ,  -j-  . . .  -J-  n* Am  ■■  2  — -  A*, 


Ax  4-  2 M,  +  ...  +  "9A.— a 


For  the  remainder  term  ^  we  obtain  the  estimate 


1***1  S 


(2»  +  2)! 

!•**  («+,)! 


(*»  +  a)! 


where,  as  before,  A  •=  h2/la2. 

The  system  found  above,  where  A  V  0,  has  a  solution,  since  its 
determinant  is  nothing  other  than  the  product  n!  times  the  Vandermonde 
determinant,  of  order  n,  composed  of  the  numbers  1,  2,...,  n,  and  there¬ 
fore  it  is  different  from  zero. 

It  would  be  Interesting  to  find  the  conditions  which  must  be  im¬ 
posed  on  A  and  at  which  the  estimates  in  Section  2  remain  valid. 

The  criteria  in  Section  2  in  the  general  case  result  in  investigating 
a  very  large  number  of  inequalities  and  therefore  will  hardly  be 
liractlcable  for  Investigating  difference  equations  with  a  large 
number  of  coefficients. 

—  4.  * 

In  the  solution  of  multidimensional  problems  new  questions  arise. 

i?Cs  -£»£  ’  "$7GP 
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We  shall  therefore  examine  below  equations  with  three  and  four 
Independent  variables. 

In  addition.  In  Section  6  we  shall  derive  difference  equations 
In  which  the  desired  values  of  u( x,  t)  In  the  nodes  of  the  segments 

t-klud  t-(k  —  (,  a, ...) 

are  related  to  each  other  by  a  system  of  linear  equations  with  a 
number  of  equations  equal  to  the  number  of  unknowns. 

Equations  of  this  type,  possessing  great  accuracy,  may  be  used 
to  determine  the  unknown  values  of  u  for  each  layer  by  solving  a  sys- 

t 

tem  of  linear  algebraic  equations,  a  circumstance  which  is  especially 
attractive,  because  it  frees  the  calculator  from  a  diverging  compu¬ 
tational  process  ( to  which  the  use  of  the  non- investigated  recursion 
formulas  sometimes  leads)  but  does  not  free  him,  of  covrse,  from 
investigating  the  error,  in  order  to  obtain  the  final  result  with 
the  desired  degree  of  accuracy. 


5.  Concluding  Remarks 


This  section  contains  a  short  critical  survey  of  certain  investi¬ 
gations  devoted  to  the  problems  of  approximating  a  one-dimensional 
heat-conduction  equation  by  means  of  difference  equations  in  the  di¬ 
rection  of  the  method  and  theory  elaborated  by  the  author  [4,  5]. 

Assuming  b  =  0  in  Relationships  ( 19)  and  (20),  we  obtain  for 
the  solution  of  the  differential  equation  describing  one-dimensional 
heat  propagation  the  finite-difference  euqation 


4 


*  *4-1  +  P  Vu »+»“  Wi*i»»  +  (a  —  a — P)  Ut,k—o, 

* 


(30) 


where: 


i 

n 


<*-2l  — 


T 


l  — 


7  — 
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the  remainder  term  has  the  form  (21).  This  same  expression  (30)  was 

,  f  ‘  ’  *  r  f  r  * :  t  * 

derived  by  us  [5i  p.  83]  using  the  method  of  undetermined  coefficients. 
Concluding  the  derivation  there,  I  remarked  that  it  is  possible  to 
make  up  a  set  of  difference  equations  of  form  (30)  for  approximating 

t 

the  heat-conduction  equation,  given  X,  then  I  wrote  out  one  such 
formula  corresponding  to  the  case  X  ■  16  and  said  that  with  its  aid 
a  final  result  may  be  attained  with  an  accuracy  of  h4,  if  the  initial 
values  of  u(x,  t)  are  known  with  the  same  or  greater  accuracy.  The 
prodf  itself  is  elementary,  therefore  it  is  not  given  there.  The  i 

j  t 

course  of  the  proof  is  the  same  as  for  the  general  linear  nonhomo- 
geneous  differential  equation  of  the  second  order  with  two  variables 
or  for  an  equation  with  three  independent  variables  (see  author's 
papers  [4,  9r  5  (Section  12)]. 

‘  ,  J 

Subsequently,  in  D.  Yu.  Panov's  manual  [7,  p.  Ill]  a  Difference 

•  (  ; 

Equation  (28)  was  derived  without  estimating  the  solution  error. 

I 

The  derivation  is  made  with  the  aid  of  arguments  similar  to  ours 

[5,  Section  12].  I 

As  is  easily  noted,  Formula  (28)  ensues  from  (30)  for  X  «  6 
and  requires  no  special  derivation, 

P.  P.  Yushkov,  studying  Formula  (30)  for  other  purposes  in 
his  early  work  [12]  also  did  not  notice  that  Formula  (28),  in  partlc- , 
ular,  was  derived  from  (30)j  but  he  subsequently  corrected  this. 

I  am  unconvinced  that  such  remarks  are  generally  desirable,  but 

t 

in  the  case  in  question  they  are  apposite  or  at  least  pardonable 
and  may  be  continued.  1  ^ j 

•  i 

- Further  study  of  difference  analogs  of  the  unidimensional  4 - ; 

heat-conduction  equation  and  the  corresponding  error  is  found,  in  1 

• - :  '  2 - •' 

particular,  in  Milne's  monograph  [8,  14]  and  in  the  works  of  hia _ i 

;  1  * 

successors  [9*  10,  11,  15]  • - 1 - 0 - ■! 

src  rnc-i  ■*  STOP  HERE 
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My  remarks  would  be  somewhat  Incomplete,  If  I  did  not  note  that 
the  proof  of  Eq.  (28)  and  Estimate  (29)  given  by  the  authors  Just 
mentioned  Is  based,  firstly,  on  a  method  of  constructing  difference 
analogs  of  differential  equations  developed  by  us  previously  [5]  and, 
secondly,  on  a  method  developed  In  other  articles  [4,  5]  for  estimating 
the  error  resulting  from  the  approximation  of  a  differential  equation 
by  a  difference  equation. 

The  estimate  obtained  by  Milne  [8,  p.  154;  14,  p.  122]  will 
completely  coincide  with  ours  (29),.  if  in  the  latter  e  and  are 

deleted,  i.e.,  no  account  is  taken,  firstly,  of  the  error  in  the 
initial  values,  and,  secondly,  of  the  rounding-off  error,  and  the 
nomenclature  is  changed. 

6.  Solving  the  Heat-Conduction  Equation  with 
an  Aid  of  a  System  of  Equations 

Up  till  now  we  have  studied  only  the  recursion  formulas  which 
allow  us  to  calculate  the  values  of  u  step  by  step.  The  present 
section  is  dedicated  to  the  study  of  the  heat-conduction  equation 
by  a  finite-*difference  method  enabling  us  to  find  the  values  of  u 
for  the  next  layer  by  solving  a  special  linear  system  of  algebraic 
equations. 

The  method  will  be  developed  for  application  to  ( 15) ,  although 
it  is  suitable  for  the  general  case. 

It  is  true  that  the  numerical  solution  of  the  heat-conduction' 
equation  by  the  method  of  thiB  section  is  considerably  more  compli¬ 
cated  than  the  solution  using  recursion  formulas,  but  on  the  other 
hand,  as  has  already  been  mentioned  at  the  end  of  Section  4,  it 
frees  the  calculator  from  having  to  analyze  the  convergence  of  the 
computational  process. 
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Consequently,  we  are  at  times  simply  compelled  to  refrain  from 
using  certain  recursion  formulas  which  make  the  computation  practically 
unsuitable  (unstable)  as  a  result  of  the  rapid  growth  of. the  error 
in  the  solution  even  when  the  error  in  the  initial  values  is  . 
imperceptible  or  as  a  result  of  rounding-off  errors. 

The  proposed  method  of  constructing  new  formulas  for  the  solu¬ 
tion  of  the  heat-conduction  equation  has  an  independent  interest, 
since,  in  the  first  place,  it  can  be  used  for  solving  multidimensional 
heat-conduction  problems,  and,  secondly,  it  requires  the  use  of 
quadrature  formulas,  in  contrast  to  methods  which  extensively  use 
formulas  of  numerical  differentiation  or  to  the  method  of  undeter¬ 
mined  coefficients  [5]. 

Among  the  many  such  formulas  let  us  dwell  on  the  following  (t 

is  considered  as  a  parameter) 

u  (*  —  A,  t)  —  iu  (x,  l)  -f  u  (*  +  A.  0=* 

«*■ —  —  h,  i)~r  •«  *'■•(*•  0-ra** (*  +  *>*)]—  (3l) 

12 

- - —  «{?(<;,  *)(*  —  h-giSx  +  h), 

240 

as  the  most  suitable  for  our  purposes;  it  was  used  by  the  author 
[16]  to  solve  eigenvalue  problems. 

Let  us  now  substitute  t  +  l  into  (31)  in  place  of  Jb,  add  the 
resulting  formula  terra  by  term  to  Formula  (3l),  and  replace  in 

_p 

the  resulting  expression  by  a  u^.  Let  us  transform  the  latter 
expression  by  means  of  a  trapezoidal  formula  of  closed  type: 

u't  (*,  t  -f  •/)  +  «* « (  v,  0  ™  ’  (  32) 

•  p 

_2_  Yu  IX,  /  +  D -- U  (r,  1)]  +  —  (x,  T,)(l<ij</r/) 

l  6 

and  of  two  more  of  the  same  resulting  from  (32),  by  replacing  x  in 
it  by  x-h  or  x  +  h.  Reducing  the  similar  terms,  we  arrive  at  the 
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difference  equation 


(6 —X)  U(r—  A,  1 4- 1)  —  (i*  t*  i<*)  U(* i  »+0  + 
+(6-X)l/(*-fA,  *+/)  +  (*  +X)U(*— A,  /)  + 
—  (r 2  —  10X)  U (x,  l)  +  ($  +  1)  U(*  +  *•  O-o. 


(3?) 


the  remainder  term  R  of  which  may  be  evaluated  with  the  aid  of  the 
inequality 


where  M«  is  the  greatest  value  of  |u^(x,  t)  J  in  rectangle  Q. 
If  now 


then  Eq.  (33)  will  be  a  trinomial  relating  the  unknown  values 

of  U  in  three  nodes  of  the  layer  t  +  Z  to  its  values  belonging  to 
layer  _t.  Similar  equations  may  be  written  for  each  node  of  the  layer 
t  +  Z.  The  equations  for  the  boundary  nodes  with  abscissas  x  *  h 
and  x  =  L-h  will  be  binomial. 

A  trinomial  linear  algebraic  system  of  equations  was  studied 
by  us  in  an  earlier  work  [17»  section  25].  There  an -analytical  method 
of  solving  a  system  was  developed,  according  to  which  the  exact 
value  of  U(h,  t)  is  found  first,  and  then,  using  the  equations  of 
the  given  system, 

U(  2h,  t)  ,  U(3h,  t)  , . . . 


are  found  by  successive  substitutions. 

The  system  derived  above  may  also  be  successfully  solved  by  using 
successive  approximations.  The  convergence  will  be  investigated 
below. 
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In  the  system  consisting  of  equations  of  type  (33)  x  and  t 
intervals  may  be‘  arbitrarily  chosen,  which  advantageously  distinguishes 
the  numerical  method  examined  from  the  preceding  one  ( Section  4)  . 

When  X  ■  6  we  again  arrive  at  recursion  Formula  ( 28) . 

It  remains  to  show  the  convergence  of  the  method  of  successive 
•  approximations  for  (33)  .  Let  us  rewrite  it  in  the  form 

[/(*,  t  +  -  A, *  +  0+  U (*  +  A.  t  +  i)J  +  A(x,  0,  / 


where 


6—X  . 

*“  12+101  ' 


(35) 


the  quantities  A(  x,  t)  depend  on  the  values  of  U  of  the  preceding  jt 
layer;  therefore  they  are  known. 

We  shall  consider  that  we  have  chosen  (completely  arbitrarily) 
the  initial  approximate  values  of  U  in  the  nodes  belonging  to  the 
segment  t  +  l .  The  result  of  substituting  these  values  into  the  right 
side  of  (34)  will,  in  general,  be  distinguished  from  u(x,  t  +  l)  ; 
we  shall  therefore  estimate  the  closeness  of  the  approximation  to 
u(  x,  t  +  l) . 

We  shall  designate  by  i(x,t  +  l)  the  difference  between 
u(x,t  +  l)  and  the  result  of  the  substitution.  Let  e  designate  the 
maximum  absolute  value  of  the  difference  between  the  exact  values 
of  u  and  its  initial  values  in  the  nodes.  We  shall  estimate  the 
difference  £(h,  t  +  l)  .  Since  u(0,  t  +  i)  is  known  to  us,  we  add  it 
to  A(x,t)  . 

Thus,  having  designated  the  error  in  the  first  approximation 
by  $(x,  t  +  l)  ,  we  shall  verify  that  at  the  moment  t  +  Z  the 
following  inequality  will  be  true  for  it 
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That  la, 

I5,(2*,<+/)1S  «(!+«)«• 

Similarly 

&(3M  +  J)1S«(i  +  «  +  «*)»• 

Continuing  these  estimates  further,  we  find  that  for  any  values 
of  k  the  following  inequality  will  be  fulfilled 

It*  (**.«+ 4!  S  —  «. 

i  — * 

if  a  <  1,  or,  according  to  (35)  ,  0  <  X  <  6  . 

Now  let  us  estimate  the  closeness  of  the  second  approximation 
to  u(x,  t  +  l)  .  As  in  the  preceding  case 

Continuing  these  estimates  futher,  we  find  that  for  any  values  of  n 
the  following  inequality  is  fulfilled 

and  the  error  £(x,  t  +  i)  will  tend  to  zero 

lim  5»(AA,  /  +0“°* 

if  a(l-a)-1,  i.e.,  if  a  <  0.5.  Since  X  belongs  to  the  interval  (0.6) 
ana  the  maximum  value  of  a  when  X  varies  in  that  interval  is  0.5, 
then  by  choosing  h  and  l  sq  that  0  <  X  -<  6,  we  ensure  the  conver¬ 
gence  of  the  successive  approximations  for  (34). 

In  order  to  illustrate  the  use  of  Eq.  (33) ,  let  us  solve 
numerically  the  differential  equation 

(36) 

dt  dx* 

0 


with  the  boundary  conditions 
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I5i(2A,  *-f /)IS«0  +«)«. 


That  is. 

Similarly 

I5X(JA.  <  +  *>IS«  (i  +  «  +  «*)*• 

Continuing  these  estimates  further,  we  find  that  for  any  .values 
of  k  the  following  inequality  will  be  fulfilled 

I6,(M.H-0IS  —  «. 

i  — • 

if  a  <  1,  or,  according  to  (35) ,  0  <  X  <  6  . 

Now  let  us  estimate  the  closeness  of  the  second  approximation 
to  u(x,  t  +  l)  .  As  in  the  preceding  case 

If,  (**,/  +  01^ 

Continuing  these  estimates  futher,  we  find  that  for  any  values  of  n 
the  following  inequality  is  fulfilled 

and  the  error  £(x,  t  +  l)  will  tend  to  zero 

lim  <+0“°* 

if  a(l-a)-1,  i.e.,  if  a<  0.5*  Since  X  belongs  to  the  interval  (0.6) 
and  the  maximum  value  of  a  when  X  varies  in  that  interval  is  0.5, 
then  by  choosing  h  and  l  so  that  0  <  X  -<  6,  we  ensure  the  conver¬ 
gence  of  the  successive  approximations  for  (3*0  • 

In  order  to  illustrate  the  use  of  Eq.  (33),  let  us  solve 
numerically  the  differential  equation 

( 36) 

dt  dbr* 


with  the  boundary  conditions 
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M«# 

and  the  Initial  conditions 


x»0)  n«o  x*l|  2 


llaCOI - * 

M 


when  t  »  0. 


Let  us  take  h  *  0.24,  l  *  0.048.  Then 

A* 

and  the  difference  equation  In  which  we  are  Interested  will  assume 
the  form 

0  (x,  t  + 1)  -  0,2  [U  (x  -  A,  t  +  i)  -f  U(x -f A.  /  +  /)]  + 

+  °,3lu(*—  A,  /)+£/(*  + A,  /)].  (  37 ) 


If  we  successively  set  In  It  x  *  h,  2h,  3h,  4h  and  t  *  0,  we  obtain 
a  system  of  equations 


^i>x“°>2  +  °«72A9a, 

0,,j= 0,2  (14,1  +  +  0,46165, 

^3,i=0»2  (y»n  +  ^*,j)  +  o,3354i, 

^«u=0»2  ^*>1  4"  ®*27^34» 


the  solution  to  which  is  given  in  Table  1.  At  the  end  of  the  table 
the  values  of  the  exact  solution  are  given  for  comparison. 

The  values  of  u  in  the  nodes  of  the  succeeding  layer  t  =  0.096 
etc.  may  be  calculated  in  similar  fashion. 


TABLE  1 


M 

ApproxlMt*  values  of  u  Valuo*  of  oxaet 
iod«»  of  lajrtr  t  »  0*04?  solution 

O  1  j 

0,92:05 

• 

0,92105 

o,U 

0,87504 

0*87597 

0*8 

o,745ii 

0745M 

•  O.72 

0,96 

0,54135 

0,2846! 

0,54158 

078462 

1.20 

0 

0 
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In  the  case  of  more  complex  boundary  conditions  the  methods 
described  above  may  be  used  almost  without  change.  It  Is  necessary 
only  to  properly  construct  the  equations,  approximating  the  boundary 
conditions. 

Thus  In  the  case  of  two  boundary  conditions  of  the  general  type 

m(o.  /)  «, ii ’,(  »,»)"  M/>.  ( 3®) 

it  is  necessary  only  to  replace  the  partial  derivatives  appearing 
In  them  according  to  the  unilateral  formulas  of  numerical  differentia 
tlon  (cf.  author's  earlier  works  [17,  18]) . 

But,  unfortunately,  when  integrating  with  the  aid  of  recursion 
formulas,  such  a  replacement  often  leads  to  considerable  distortion 
of  the  solution,  since  In  the  actual  calculation  of  the  partial 
derivatives  the  already  calculated  values  of  u(x,  t) ,  which  were 
subject  to  error,  ax*e  used. 

Therefore  It  Is  desirable  to  find  some  new  means  for  approxi¬ 
mating  the  boundary  conditions  when  solving  equations  with  the  aid  of 

m. 

It  is  natural  to  expect  that  they  may  be  obtained  by  using  for 
the  approximation  both  the  given  initial  conditions  and  the  given 
differential  equations  simultaneously.  Let  us  now  turn  our  attention 
to  the  construction  of  these  formulas. 

We  shall  begin  with  a  consideration  of  the  quadrature  formula 

»(*—*,  0  —  « (*  —  A,  l)  =  ihu' t  (x  —  k,  t)  +  ^  jy) 

+Y  A’««  (r  ^  h,t)+  A*uirt  (*,  t)  -h  R, 

where 


I*! 


2h* 
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Here  Ms  denotes  the  maximum  value  of  !*£(*»  01  in  the  region  Q.  A 
derivation  of  this  formula  appears  on  page  1210  in  one  of  the  author's 
previous  articles  [19]. 

The  method  of  constructing  Eq  (33)  with  the  aid  of  (31)  can 
be  applied  almost  without  change  to  Eq.  (39) •  After  employing- 
trapezoidal  formula  (32),  it  is  only  necessary  to  replace  x  by  h» 
in  order  to  obtain  the  final  formula.  We  thus  obtain  the  relation¬ 


ship: 

+  ±  a  )  » (o,  /  +  i>-  -  aA  l**. (o.  <  +  0  +  «•’. (o.  01 + « (»*.** + 0  + 
|| (3a,  /)  -  l«  (A.  /  +  /)-« <*.  01-  ( *  *>  r» 


(40) 


the  remainder  term  R*  of  which  may  be  evaluated  with  the  aid  of  the 
inequality 

I#*  |  •  ‘ 

45  > 


where  M3  is  the  maximum  value  of  /)  in  the  rectangle  Q,  and 


Determining  the  partial  derivatives  u',(o,  u',(n,  i)  for 

a i  /  0  from  the  upper  equation  of  (38)  and  substituting  them  into 
(40)  leads  to  the  equation  in  which  we  are  Interested. 

It  is  also  possible  to  derive  exactly  an  equation  similar  to 
(40)  containing  nVl,  /+/>□*  v'.{L,n.  and  then  eliminate  from 
it  the  derivatives  with  respect  to  the  spatial  variable,  employing 
for  this  purpose  the  lower  equation  in  (38)  .  It  should  be  noted 
only  that  before  we  proceed  to  construct  the  formula  in  which  we 
are  interested  we  must  replace  h  by  -h  in  Quadrature  Formula  (39) 
and  then  set  x  =  L  -  h  in  the  formula  thus  obtained. 

As  an  example  we  shall  again  solve  Differential  Equation  (36) 
with  the  initial  condition 
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"  (x •  o);  30  COR  —  x  (o  t't  x  I  )t 
•  2 

and  the  boundary  conditions 

«’•(<>.  /)=o,  «(i,  /)-o, 


,  i.e.,  we  shall  solve  the  problem  of  the  cooling  of  a  rod  with  a 
thermally  insulated  lateral  surface,  when  the  Initial  temperature 
distribution  is  known,  assuming  that  one  of  Its  ends  is  thermally 
Insulated,  while  the  other  is  kept  at  a  constant  temperature  0. 

Let  vs  take  a  =  1,  h  =  0.1,  and  1=1=  1/120;  we  then  obtain: 
X  =  1.2  and  the  difference  equation  will  again  be  of  type  (37). 

By  setting  in  it  x  =  h,  2h,  . . . ,  9h,  successively,  and  t  »  0,  we 
obtain  a  system  of  nine  equations  in  ten  unknowns  u.  (k  *  0,  1,..., 
9)* 

yj,«  =  o.2(r,„  +  11,76634, 

4.  Lji)  +  *1. 27217, 

,0i5<’°44» 

=  9,58868, 


the  last  of  which  satisfies  the  boundary  condition  u  (l*y|§)  »  0. 

We  obtain  the  missing  equation  after  satisfying  the  boundary 

condition  in  ( 40)  u^(0,  t)  =  0  at  the  points  (0,  0)  and  (0,  * 

Rejecting  the  remainder  term  and  taking  h  =  0.1,  we  obtain  the  follow- 

P».i  “*°,s  (Uitl  -j-  Ultl)  +  8,38082, 

ing  formula  ,.  ...  .... 

^m*”0*2  +  U  t,i)  +  6,96659, 

^••i“°»2  (^»ii  4* 

_  +  1.85410, 

which  is  accurate  up  to  h5.  * 

We  shall  give  the  results  of  the  computations  (Table  2).  The 
values  of  the  exact  solution  are  given  at  the  end  of  the  table 
for  comparison. 
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TABLE  2 


Solution  of  Eq.  (?6)  for  nodes  of  the  layer 
t  -  1/120 


Initial  *alu»»  of  » 

Approx  inoti  vmluta  In 

vauoa  If  tho 

in  tho  nodoo  of 

nodta  of  tho 

X 

tho  loor* r  t  •  0  . 

t  •  1/130 

ixaot  solution 

o 

20 

19,5930! 

19.59300 

0,1 

*9.75377 

19.35*74 

>9.35*78 

0,2 

19,02113 

18,65401 

18,63405 

*3 

17.82013 

1  *7.45745 

1M5749 

0,4 

16,18034 

1585103 

15.8;  107 

0,5 

o,6 

14,142*4 

m.75571 

13,8543* 

11,51646 

13.85435 

11,51648 

o,7 

■  9,07981 

8,89502 

8.89504 

0,8 

6,18034  ] 

6,05458 

6,05457 

o,9 

3.1*869  j 

3.06501 

3,06502 

1,0  ■ 

0  1 

0 

0  ' 

The  system  was  solved  by  the  iteration  method. 

The  convergence  of  the  iteration  can  be  proven  for  the  general 
case,  i.e.,  for  a  system  of  equations  consisting  of  i/h  -  1  equations 
of  type  (3*0  and  one  more  equation 

17  (o  /+/)=» _ ? - f  £7(sA, /  +  /)  —  — IT  (A, /-{-/>  0. 

3+4*L  3  3 


which  results  from  ( 40) .  Here  B(h,  X,  t)  is  known,  since  it  is  a 
function  of  the  values  of  u  in  the  preceding  layer. 

The  following  inequality  will  hold  true  for  the  error  in  the 
first  approximation  at  the  moment  t  +  l 


i-«  3  +  4I 

provided  that  a  <  1  and  consequently  that  0  <  X  <  6,  in  accordance 
with  ( 35) .  Continuing  on,  we  find  that 


Thus  the  convergence  of  the  interations  is  ensured  by  the  condition: 
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From  the  amount  by  which  the  left  side  of  this  inequality  differs 
from  unity  it  is  possible  to  Judge  the  rapidity  with  which  ln. 

When  X  *  1.2,  the  left  side  of  (4l)  is  equal  to  is  21/52,  and 
consequently  there  is  convergence.  This  value  of  X  is  not  the  best 
of  the  number  of  possible  values  ensuring  rapid  convergence  of  the 
iterations  in  our  problem. 

In  order  to  approximate  the  boundary  condition  u^(0,t)  «  0 
it  is  also  possible  to  use  Taylor's  formula,  the  given  differential 
equation,  and  Formula  (52) .  We  obtain! 

7T  _  ^t>»  +  k\u  -j-  (X  —  l)  17.* 

U*1  - 7+1 - * 

Then  for  the  example  considered  above  we  arrive  at  the  formula 

V  +  UyJ 

1,1  aa  ’’ 

which  is  accurate  up  to  h3. 

7-  Multidimensional  Problems 

In  solving  multidimensional  problems  questions  arise  which  re¬ 
quire  explanations.  Further  examinations  must  be  carried,  out 
separately  for  two-dimensional  and  three-dimensional  space. 

We  shall  first  examine  the  propagation  of  heat  in  a  thin  plate 
with  a  contour  of  arbitrary  shape  y.  For  this  case  we  have  the 
differential  equation 


where  A  denotes  the  Laplacian  operator! 

*  # 

is*r+*- 

u  is  the  temperature  of  a  certain  point  (x,  y)  at  the  moment  t,  and 
a*  is  the  thermal-diffusivity  coefficient. 

We  are  thus  confronted  with  the  problem  of  the  propagation  of 
a  temperature  u(x,  y,  t)  in  a  finite  principal  region  D  bounded  by 
the  curve  y  .  The  problem  consists  in  finding  a  solution  to  Eq.  (42) 
which  satisfies  the  initial  condition  and  the  boundary  conditions. 

Let  us  first  examine  the  boundary-value  problem,  when  the 
initial  temperature  distribution  in  the  plate  is  known 

»  (*,  y,  <*)  -  ?  (x,  v)  (  ^3) 


and  during  the  whole  time  of  observation  the  temperature  on  the  edge 
of  the  plate  y  is  kept  equal  to 

i* (s,  r,  y,  /),  (  44) 


where  f(x,  y,  t)  is  a  given  function  on  the  curve  y,  and  r  and  £  are 
the  coordinates  of  a  variable  point  on  this  curve. 

Equations  (42)  with  Conditions  (4j5)  and  (44)  has  been  examined 
previously  [4,  5].  The  author  [4]  has  examined  an  algebraic  analog 
which  approximates  the  equation  vJ  ih  an  accuracy  up  to  h2,  where  h 
is  a  side  of  a  square  of  the  network,  and  in  another  paper  [5]  an 
analog  which  approximates  (42)  with  an  accuracy  up  to  h4  was  given. 
The  latter  analog,  i.e.,  the  improved  difference  equation  ([5]  P-  92) 
approximating  ( 42) ,  has  the  form! 

U(x,  +  U{x+h,  y+h,  1)+  U(x  —  h,  v  +  M+  (45) 

36  l 

+  U(x-h,  y-h,  t)-i-r(x  +  h,  y-h,  Q  +  4lP(x  +  A,.*0  + 

•+  V  (x.  y~h,  <)  -f  i:{x  -  h,  y,  I )  -i  V  (x,  >•  -  h,  <)]  +  1 61’  (x,  y,  t)  1. 
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where  x  and  2.  are  the  coordinates  of  the  main  node  of  the  square 
network,  i.e.,  that  (fixed)  node  (x,  y)  for  which  Equality  (45)  la 
written,  and  i  is  the  spacing  with  respect  to  t. 

For  convenience  Formula  (45)  is  presented  In  the  form: 

4  * 


(46) 


where  2  (respectively  2  )  is  the  sum  of  the  values  of  U  at 
k  =  1  k  =  5 

the  moment  t^  in  the  network  nodes  at  a  distance  h  ( respectively 
*V2~ h)  from  the  main  node  0 . 

Considerably  later.  Formula  (45)  was  derived  by  an  operator 
method  by  W.  E.  Milne  ([8],  p.  150;  [14],  p.  157)  and  written  in 
symbolic  form  with  the  aid  of  a  pattern 


1 

;  1 

4 

x  . 

U(x,y,t  +  l)-*±.  4 

16 

4 

30  j 

■  ,  UL 

4 

1 

U  (*i  y,  <)•* 


(47) 


Then  Eq.  (47)  was  proven  again  ([11],  p.  152;  [15],  P*  118)  and 
the  proof  was  completed  with  the  remark  that  it  could  be  conveniently 
used  for  solving  the  heat-conduction  equation  in  two  spatial  variables 
on  computing  machines. 

Hexagonal  networks  (these  are  known  as  triangular  networks 
in  the  literature)  for  approximating  Eq.  (42)  have  been  studied 
by  P.  P.  Yushkov  [20]. 

In  none  of  these  papers  do  we  find  an  estimate  of  the  error  in 
the  numerical  solution.  Therefore  Eq.  (46)  is  generalized  in  this 
section  for  rhombic  networks,  and  an  estimate  of  the  error  in  the  nur 
merical  solution  of  the  heat-conduction  equation  found  with  the  aid 
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of  this  generalized  formula  Is  given.  The  generalized  formula  leads. 
In  particular,  both  to  Formula  (46)  mentioned  above  and,  consequently, 
to  the  Identical  Formula  (47) ,  and  also  to  a  countless  number  of  other 
formulas  which  are  of  practical  Importance  In  certain  cases. 

All  of  these  formulas  are  Important  In  studying  the  propagation 
of  heat  In  square,  rhombic,  and  hexagonal  plates,  although  they  may 
also  be  successfully  used  for  plates  with  any  boundary  configuration. 
In  these  cases  it  Is  only  necessary  to  add  to  formulas  of  type  (46) 
special  formulas  [6]  for  points  adjoining  the  boundary.  Afterwards, 
the  method  developed  by  the  author  [6]  for  curvilinear  boundaries 
was  used  in  several  books  [8,  14],  §  65;  [11],  15],  §  8,  6;  [21], 

§  12. 


; 


A.  Let  us  now  proceed  to  derive  a  formula  for  approximating 
(42)  for  the  case  of  rhombic  networks.  We  shall  base  the  derivation 
on  a  previously  developed  method  [22],  §  which  involves  constructing 
the  relationships  between  the  values  of  the  functions  of  two 
variables.  This  will  permit  us  to  write  out  the  expansion  of 
u(x,  y,  t)  for  the  main  node  (x,  y) ,  which  lies  in  the  center  of  a 
network  rhombus  with  sides  2h  (Fig.  l)  and  diagonals  parallel  to  the 
coordinate  axes.  The  expansion  has  the  form! 
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(48) 


«*t*  (*>•?»  0“ 


QiUi + Axh*  l  v  (*.  y,  0  +  A*  k  (*.  y>  0  + 


where 

t 

— 2sin,d>608,*>» 

14  +  -f  3  ctg* «, 

8 

£  eA-4(«,+N+«.+'0+(s*i,»-*)-2“+ 

*=*I 

-f  (3o»tr*u>— i) 

2 

|/J|  <  [j j  sin*  «o  cos’  w  —  4  (sin*  u>  -f*  cos*  w)  +  (sin  «  +  cos«)*JW«» 
4S 


(49) 

(50) 


where  Uj,  U2,  U3,  and  U4  denote  the  values  of  u(x,  y,  t)  at  the 

moment  at  the  vertices  of  the  rhombus  [with  center  In  the  main  node 

(x,  y)  ]  with  sides  2h  and  diagonals  parallel  to  the  coordinate  axes 
( the  odd  subscripts  refer  to  vertices  lying  on  the  horizontal  diagonal* 
and  the  even,  to  those  lying  on  the  vertical  diagonal) \  u5,  u6,  uT, 
and  Us  denote  the  values  of  u(x,  y,  t)  at  the  mid-points  of  the  sides 

of  the  same  rhombus tc  is  the  value  of  the  angles  formed  by  the  hori¬ 

zontal  diagonal  of  the  rhombus  with  the  sides  of  the  rhombus;  and 
Me  denotes  (here  and  henceforth)  the  upper  bound  of  the  values  of 


6%u 

d*u 

d*u 

d*u 

dx* 

> 

dx'df 

|  dx*  dyk 

d? 

within  a  certain  three-dimensional  region  of  space  x,  £,  and  jt  in 
which  the  solution  is  sought. 

In  order  for  the  coefficients  of  to  be  non-negative,  in  all 
that  follows  ad  will  be  determined  by  the  condition! 
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Prom  (49)  it  follows  that 


o 


(52) 


Then  if,  together  with  ( 48)  and  ( 42) ,  we  make  use  of  the  Taylor 
series  expansions  of  u(x,  y,  t  +  l)  and  u(x,  y,  t  +  2l)  exactly  as 
we  did  in  the . one-diemensional  case  (§4),  we  obtain! 

*  “  (*.  y,  t -f 1)  -f  (x,  y,  t  +  2t)  -f  (a,  —  a  —  0)  u (*,  y,  /)« 

8 

—  «<«<  4-  R  +  «  /?,  +  ?RV 

1=1 

where 


(53) 


2j1tX*  —  l2X  —  4X*8iDt«i>C0StM, 

?=o,5  —  jl  -j-  2X*  siii* <■> cos*  <0, 

,•  A* 


(54) 


Uf  ’ 

while  for  Rj  and  R2  we  have  the  inequalities 


(55) 


3*’ 


Ji*.  * 


Let  us  now  choose  X  such  as  to  make  0  equal  to  zero.  This  will 

allow  us  to  discard  the  term  0u(x,  y,  t  +  21)  in  (53),  and  we  shall 

« 

arrive  at  a  recursion  formula  allowing  us  to  determine  the  values  of 
u  at  the  moment  t  +  l  from  its  values  at  the  moment  ^t. 

Moreover,  we  find  from  Eqs.  (54)  that 


Sin*  CD  cos’  <D 


X«. 


i»5 


sin”  cd  cos”  cd 


*  a  —  fi  •  8  - 


sin*  CD  cos1  CD 


Consequently,  when  P  =  0 


„  z  .  j  k  6—8  sin*  cd  cos*  cd 

u  (*>  + 0  - - » (*,  y,  0  + 


.  s 

sin*  cd  cos*  cd 


\  ^  ,  D  ,  sin*  cd  cos*  m 

/  <  «c«c“r/?»  + - *• 

i*i  * 


It  is  easy  to  see  that 


(56) 


C — Ssin'wcoj*  «><>. 
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We  also  know  that  all  the  a^'s.in  (49)  are  greater  than  or 
equal  to  zero.  If  to  lies  In  ( 51)  .  Consequently  variation  of  to  be¬ 
tween  ^  and  ^  entails  the  non- negativity  of  the  coefficients  of 
Eq.  (56).  In  all  that  follows,  therefore,  let  us  agree  to  consider 
Formula  ( 5 6)  for  those  values  of  ®  for  which  Inequality  ( 51)  Is  valid. 

Thus  In  order  for  the  computational  process  generated  by  Fbrmula 
(56)  to  converge.  It  is  enought  that  the  relationship  between  the 
h  and  l  Intervals  be  taken  such  that 

In  any  case,  if  ®  belongs  to  (51)  ,  we  obtain  by  using  (56)  and 


(52) 


2  sin* w cos* «  , 

— + - ; - K  — 8)-l, 


At, *+>=  i>  8=1, 


and  consequently  the  error  In  the  numerical  solution  to  ( 42)  found 
with  the  aid  of  (56)  (without  the  remainder  term)  In  any  node  of  the 
rhombic  network  at  any  moment  t^  T  will  satisfy  the  Inequality 


where 


™ W *  H|. 


where  Inequalities  (50)  and  (55)  should  be  fulfilled  for  Rx  and  R. 

Now  let  the  boundary  contour  7  be  the  network  contour,  i.e.,  it 
consists  of  the  nodes  and  links  of  the  rhombic  network.  Then  in  the 
estimate  Just  obtained  e  should  be  tkan  to  mean  the  absolute  value  of 
the  error  in  the  function  u(x,  y,  t)  (0^;  t^  T)  in  the  nodes  of  7. 
Recalling  the  estimates  for  R!  and  R,  we  obtain  the  following  theorem! 
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If  the  solution  of  Eg.  (42)  with  Conditions  (43)  and  (44)  Is 


sought  with  the  aid  of  Formula  (56)  (without  the  remainder  term) , 
then  the  error  In  any  node  of  the  network  ( lying  Inside  7)  at  the 
moment  of  t  +  T  satisfies  the  Inequality. 


|;(r,  v.  T)  1-5  e 


Ca*T\lth* 
270  ’ 


(57) 


where 

G*»  K'O  Mil4  (•>  COS*  (1)  -f  l2SiU*«i>C08*M  4(piu‘w  +  cos'**)-}- 

\ 

-f  (sin  w -f- cos w)‘ 

In  the  case  of  a  square  network  with  the  diagonals  of  the  squares 
parallel  to  the  coordinates  axes  o  =  i.e.,  X  ■  6,  from  it  we  can 

deduce  the  following  theorem  concerning  the  estimate  of  the  error  in 
the  solution! 

Theorem.  Let  the  solution  to  Eq.  (42)  with  Conditions  ( 4  3) 
and  ( 44)  be  sought  with  the  aid  of  the  formula 


V(x,  <  +  0=-~  lt7i,»  +  t'„i  +  tr4,i+ 

3* 

_h  •+■  Ut,t  t  V iit  4-  «)  -f-  16  Ut,  *1, 


(58) 


where  UQ>^  is  the  value  of  U(x,  y,  t)  at  the  moment  t^  at  the  center 
of  the  network  square. 


Then  the  error  in  any  network  node  ( lying  inside  7)  at  the 
moment  t  =  T,  in  accordance  with  (57) ,  satisfies  the  inequality 


saTAf,*4 

27 


(59) 


where  h  is  the  side  of  a  square  (inclined  at  an  angle  of  45°  to  the 
x-axls)  . 

The  theorem  is  proven  in  the  same  way  for  hexagonal  networks 
with  a  network  contour  7  consisting  of  the  nodes  and  links  of  the 
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network.  For  such  a  network  the  following  theorem  applies. 

Theorem.  If  the  solution  to  Eq.  (42)  with  Conditions  (4?)  and 
44)  is  sought  with  the  aid  of  the  formula 


n  6l/t,  i  -f-  -f*  Uj,  *  -M4 .« + U«.  •  +  Uw  +  U$, « 

iW - r - u - - - ’ 

where  Uc  t  is  the  value  of  U(x,  y,  t)  at  the  moment  Jfc  at  the  center 


of  a  certain  hexagon  of  the  network,  while  a  k  *  : 

8)  is  the  value  at  its  vertices,  then  the  error  in  any  network  node 


( lying  inside  7)  at  the  moment  t  =  T  satisfies  the  inequality 

li(., ,,.)!<.+  ‘"'.Iff-.  (6°) 


where  h  is  the  side  of  a  hexagon  of  the  network. 


prove  this  theorem,  it  is  enough  to  take  <0  *  ^  in  (56)  and 


(57)  • 


For  the  sake  of  completeness,  let  us  consider  a  square  network 
with  the  sides  of  the  squares  parallel  to  the  coordinate  axes.  Let 
us  return  to  Formula  (48)  and  in  it  set 


8  _4, 

•  f-  X/.*.-4E*>  +  E"- 


*r=I  *  =  S 


in  accordance  with  Formula  ( 13)  derived  in  chapter  I  ( §  2)  in  our 
earlier  publication  [5].  For  this  case  Estimate  ( 50)  is  already 
invalid.  Calculation  shows  that  it  should  be  replaced  by! 

uiaizM.. 

4S 

Let  us  now  examine  the  formula  thus  obtained,  together  with  (53)  and 
the  left-hand  equalities  in  (54) ,  when 


2 

We  then  obtain  X  *  6  and  a  *  36,  so  that 
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For  these  values  of  X,  a,  and  M  we  arrive  at  Formula  ( 46) ,  for  which 

+  .  (61) 

10  • 

where  h  is  the  side  of  a  square  of  the  network. 

As  is  apparent.  Formula  (46)  and  (48)  have  approximately  the 
same  degrees  of  accuracy. 

Thus  in  view  of  Estimates  ( 59) ,  (60),  and  (6l)  we  arrive  at  an 
interesting  result,  namely,  that  out  of  the  three  formulas  with 
identical  h’s  Just  examined  [for  the  numerical  solution  of  Eq.  (42)] 
the  formula  for  hexagonal  networks  generally  yields  the  most  accurate 
result. 

Let  us  now  proceed  to  a  consideration  of  the  general  case.  Let 
us  determine  a  and  £  in  (54)  with  the  requirements  that  for  any 
rhombic  networks  characterized  by  the  two  numbers  h  and  cd  and  the 
time  interval  i,  l.e.,  for  any  X  which  figures  in  (54),  the  in¬ 
equalities 

—  o,  -«,  +  *  +  ?—  o* 

are  fulfilled.  In  accordance  with  the  discussions  in  §  2,  this  is 
necessary  for  the  convergence  of  the  computational  process  and  for  the 
estimate  of  the  accuracy  ahcieved. 

Using  ( 54) ,  we  can  verify  that  the  first  and  last  of  these  three 
inequalities  are  fulfilled  simultaneously  when 

X‘> - * - 

ftio'ttcos1* 

for  all  values  of  o>  In  the  Interval  (5l)j  moreover,  a  study  of  the 
second  Inequality  indicates  that  it  can  be  realized  for  any  co  in  (51) 


when 


9— t 


4  sin* »  cos*  w 


3  *  S 


•-9  +  t  * 

4  sin*  u>cos*w 


where 

T"K  57  —  . 

Furthermore,  no  matter  what  the  value  of  <u  is  from  (51),  the 
last  Inequality  may  be  replaced  by  the  Inequality 


3  _ 9  +  ~  (62) 

sin*  w  cos* »  4  sin*  w  cos  w* 

Thus,  X  and  cu,  determined  by  Inequalities  (51)  and  (62),  satis¬ 
fy  all  three  Inequalities  considered  above  and,  consequently, 
guarantee  the  convergence  of  the  computational  process  carried  out 
for  a  rhombic  network  with  the  aid  of  Formula  ( 55) ,  In  which  It  Is 
necessary  to  take  a  and  P  In  accordance  with  Formulas  .(54)  »  to  re¬ 
place  the  sum  Z  a.,u.,  in  accordance  with  (49),  to  substitute 
1  1  1 

,4  -I-  5  ts*  w  +  }  cot8  <0-91  +  21*  sin* » cos*  <0, 

In  place  of  the  quantity  aQ  —  a  —  p,  and  to  solve  the  equation  ob¬ 
tained  for  u(x,  y,  t  +  2l). 

It  remains  to  give  estimate  of  the  error  In  the  numerical 
solution  sought  with  the  aid  of  the  recursion  formula  Just  obtained. 
For  this  purpose,  we  should  evaluate 

M=rill*:+*l*,l+?l*,l).  (65) 

using  Inequalities  ( 50)  and  ( 55) .  From  this  there  results  the 
following  theorem! 

If  the  solution  to  Eq.  (42)  with  Conditions  (43)  and  (44)  Is 
sought  with  the  aid  of  recursion  Formula  ( 53)  solved  relative  to 
u(x,  y,  t  +  21)  ,  then  the  error  In  any  node  of  a  rhombic  network 
(lying  Inside  7)  at  the  moment  t  «  T  satisfies  the  Inequality 
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I?<*.  r,t)\%t  +  L(M  +  V), 

♦  * 

provided  that  a  and  P,  figuring  In  (63) t  are  determined  by  Eqs.  (54) , 
In  which  X  and  cd  must  be  taken  In  accordance  with  Inequalities  (51) 
and  ( 62)  . 

Here  c  also  denotes  the  maximum  absolute  error  with  which 'the 
values  of  u(x,  y,  t)  are  calculated  In  the  nodes  of  the  network 
contour  7  when  t  ^  T  and  In  Its  Internal  nodes  when  t  *  0  and 
t  -  l. 

If  these  values  are  computed  to  an  accuracy  of  h4,  then  we  obtain 
the  values  of  U  at  the  moment  T  with  the  same  accuracy. 

Setting  a)  =  we  obtain  from  (62)  for  a  square  network  (with 
the  diagonals  of  the  squares  parallel  to  the  coordinate  axes)  the 
Inequality 

12SXS5-H/47  (64) 


(65) 


and  the  computational  process  carried  out  with  the  aid  of  the  formula 

t/„  i+,,-(o,j  X*  -  >X)-1  ((X*  -  1 2X)  {/„ ,+,  -  (o,s  X»  — 9X  +  20)  + 

4  « 

Uk,d,  • 

which  results  from  (53)  >  will  converge  if  X  Is  contained  In  the  Inter¬ 
val  of  (64)  . 

As  soon  as  X  satisfies  (64),  Eq.  (65)  leads  to  values  of  TJ 
subject  to  errors  determined  by  the  Inequality 

15  (*.  y,  01 S*  +  —  •*  TM#. 

9  0,5  X*— 3  X* 

If  we  now  set  X  *  12,  we*  may  write  Eq.  (65)  in  the  form 

4  • 


3b 
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In  this  case  the  error  £(x,  y,  T)  in  the  internal  nodes  of  the 
network  contour  at  the  moment  T  satisfies  the  Inequality 

I5(*.y.  T) + 

27 

t 

The  method  set  forth  above  for  solving  Eq.  (42)  is  also  appli¬ 
cable  to  curvilinear  boundaries.  In  this  case  it  is  only  necessary 
to  set  up  and  use  special  formulas  for  the  points  adjoining  the 
boundary,  exactly  as  was  done  in  our  previous  publication  [6,  §  7]» 

B.  Using  Formulas  (7)  and  (8)  from  the  author's  earlier  paper 
[23],  we  may  in  a  way  similar  to  what  was  done  for  the  one-dimen¬ 
sional  case,  derive  a  system  of  five-term  linear  algebraic  equations 
for  the  numerical  solution  of  (42).  The  solution  of  such  systems 
requires  fairly  tedious  work,  if  h  is  small,  and  the  number  of  equa¬ 
tions  is  consequently  large,  even  If  the  method  of  successive  approxi¬ 
mations  is  used.  However,  the  possibility  of  using  the  method  for 
any  spacing  relationships  makes  it  nonetheless  practical.  We  shall 
therefore  describe  the  method  in  its  general  outlines. 

Let  us  rewrite  Formula  (7)  of  the  aforementioned  paper  [23] 
in  the  form 


20  |“4 


48  4^ 

y>  ^  ^  tUkji  —  °»5  ^*r  ^  <  “f“  E  a"*’4  (66) 

ksM  1  4  =  ^  4*1  ^  • 


where  u^.  t  denotes  the  values  of  u(x,  y,  t)  in  the  specially  selected 
nodes  of  the  square  network  [used  for  ( 46)  ]  at  the  moment  t^.  Formula 
(66)  is  accurate  up  to  he. 

Substituting  t  +  l  into  (66)  in  place  of  t,  we  add  the  formula 
obtained  to  Formula  (66),  replace 


1,  a,  3.  4) 
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by  a"*a't(Xjc,  y^,  t)  in  the  sum,  and  use  Eq.  (52),  which  is  written 
for  two  spatial  coordinates. 

Thus  we  obtain  the  equation 


(20  +  8X)  (/,,  i+i  —  (4  -  -  X)  V  14,  i+i—  /  14.  i+i  * 

±-s 

4  « 

—  —  (20  8X)  Ut,  «  -f-  (4  +  *)  /  14. « -}-  »• 


(67) 


After  writing  out  these  equations  for  all  the  internal  nodes 
of  a  square  plate,  we  arrive  at  a  system  consisting  of  n  linear 
equations  in  n  unknown  values  of  the  function  u( x,  y,  t  +  i) ,  where 
jn  is  the  number  of  internal  nodes  of  the  square  bounding  the  plate. 

The  system  thus  obtained  may  be  solved  by  the  iteration  method. 
For  the  proof  of  the  convergence  of  the  iteration  process  for  any 
initial  approximations  of  t  +  f  it  is  sufficient  to  require  that 


In  practical  computations  we  shall  assume  that  the  values  of  u. 


k,  t 


already  computed  on  the  preceding  layer  are  taken  as  the  initial 
approximation  for  t  + 

We  shall  rewrite  (67)  in  the  form 


where 


4  " 


1 


5 


(68) 


•4 


•  20  -p  8X  20  8X 

48 

At*=b[( 4  +  X)  £]  17*.  1-  (ao  —  8/.)  UmI- 
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Thus,  the  quantities  A^.  are  linear  combinations  of  the  values 
of  U  at  the  nodai  points  at  the  moment  t,  combinations  which  shall 
hereafter  be  assumed  to  have  been  computed,  l.e.,  the  values  of 
in  (68)  will  hereafter  be  assured  to  be  known. 

We  shall  now  take  the  numbers  t  as  the  initial  ( zero) . 
approximation  and  set 

lcr«.t+<-  -  tVi-> 

throughout  the  entire  square  region.  The  values  of  are  known 

to  us  in  the  boundary  nodes,  and  therefore  the  difference  between 
Uk  t  +  i  an<*  It®  zero  approximations  will  be  assumed  equal  to  zero 
in  the  boundary  nodes. 

Let  us  now  compute  the  values  of  U  in  the  network  nodes  closest 
to  the  boundary  of  the  square  7,  using  for  this  purpose  the  values  of 
u  given  on  7  at  the  moment  t  +  l  and  the  initial  system  of  values 
of  for  all  the  other  nodes  encountered  in  (68) .  In  this 

way,  at  least  three  values  of  U  in  nodes  lying  on  7  will  figure 
In  each  of  Eqs.  (68).  Therefore,  after  denoting  the  error  of  the 
first  approximation  by  £i(x,  y,  t)  ,  we  shall  verify  that  at  the  moment 
t  it  satisfies  the  inequality 


!;,(*.?.  /  +  0.S(}«  +  2*)*=-«»*  (69) 

where 

2'*  ■+■  8X  * 

Thus  Estimate  (69)  holds  true  along  the  square  7,,  on  which  the 
boundary  nodes  lie. 

Then,  using  the  values  obtained  for  t+l  in  the  nodes  of  yt, 
the  initial  system  of  values,  and  Formula  (68),  we  calculate  the 
values  of  U  in  the  nodes  of  the  boundary  72  of  the  square  which  is 
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next  in  proximity  after  7j.  For  the  error  in  the  first  approximation 
at  the  nodes  of  y2»  the  estimate 

r2*)«  (3a  +  aJ)»  - " a(i 4-*)*. 

will  be  valid.  Analogously 

!?»(*•  .r.<  +  01  '-c  («  +  a*)«(i  +*)«  +  (J«+2*)*S 

:'«(i  + 

where  4i  now  denotes  the  error  in  the  first  approximation  in  the  nodes 
of  the  boundary  y3  of  the  square  which  is  next  in  proximity. 

Continuing  these  estimates  for  (68)  with  the  non-negative 
coefficients  a  and  b  until  all  the  nodes  of  the  main  square  are 
exhausted,  we  verify  that  the  estimate 


I  —  a 

finally  holds  for  the  first  approximation,  since  a  <  1‘ for  any  X^  4: 

Using  the  first  approximations,  we  calculate  the  second  approxi¬ 
mations,  etc.  by  means  of  Formula  (68)  in  a  way  similar  to  that  used 
when  we  calculated  the  first  approximations,  etc.  The  error 
estimate  for  the  m-th  approximation  has  the  form: 

|U<*.j.!+0ls(ri;)V 

whence  follows  the  convergence  of  the  successive  approximations  when 

'7 

The  most  rapid  convergence  occurs  when  X  *  4.  In  this  case 
Formula  (67)  assumes  the  simple  form: 

8 

$a  '  l/». i+4 -f  8 

i  -s 


£*•+£ 


(70) 


Using  ( 70) ,  we  arrive  finally  at  a  system  consisting  of  five- 
term  equations  relating  the  unknown  values  of  u  in  five  network  nodes 


at  the  moment  t  +  l  to  Its  nine  known  values  at  the  nodes  at  the  mo-  i 
ment  t. 

We  obtain  another  formula  for  (42)  by  using  Formula  (8)  from 

\ 

the  author's  aforementioned  paper  [23].  Calculation  yields  the  formula!  i 


(20  +  iol) 


f»iM  + 


Ut,  1+1  -f-  (l — o,$l) 

+4/”*  + 0,5!) y  Ut,  1 — (20 — 10 x) c/»t <» 

4=1  irt 


which  assumes  the  simple  form 


4  4  8 

40  U$tt+,w* 4  ^  tV«+i  +  4  2^  2^ 


for  X  =  2. 

These  formulas  are  accurate  up  to  h®.  Finally  It  Is  possible 
to  prove,  although  we  shall  not  take  the  time  to  do  so,  however, 
that  If  the  general  formula  of  the  preceding  paragraph  Is  solved  for 
uo,  t+l  and  an  Iteration  Is  formed  with  the  aid  of  the  formula 
obtained,  then  convergence  will  take  place  for  all  values  of  X  for 
which 

j-<X§!4. 

The  method  may  be  applied  to  regions  with  curvilinear  boundaries; 
in  these  cases  it  is  only  necessary  to  set  up  and  use  special  formulas 
at  the  points  adjoining  the  boundary. 

Thus,  using  Formula  (6.2)  in  the  literature  [22]  for  the  layers 
t  and  t  +  l ,  Eq .  ( 42)  ,  and  a  formula  of  type  ( 32)  written  for  two 
spatial  coordinates,  we  obtain  for  the  boundary  points  the  formula 
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/  »  ,  1  ,  '•  \r  1  (  £}>  '44  ,  UvtU  \ 

V  M»~  ‘  V,  '  A,  +  *,l  A,'  +  A,  J 


.where  n  is  the  value  of  U  at  the  moment  t  at  a  boundary  point  with 
the  number  v,  while  Uj,t,  U2,t,  U3,t,  U4,t  are  the  values  of  U  in 
the  nodes  separated  from  the  main  node  ( the  node  with  the  number  v) 
by  a  distance  of  hi  (to  the  left)  and  h2  (to  the  right)  in  the  direc¬ 
tion  of  the  horizontal  and  h3  ( upwards)  and  h4  ( downwards)  in  the 
direction  of  the  vertical,  respectively;  the  values  of  the.  function 

u(x,  y,  t  +  l)  ,  denoted  by  U  t+l  and  II  t+J(k  “  1»  2*  *0 

h2 

have  the  same  meaning.  As  for  X,  it  is  equal  as  before  to  -  , 

la2 

where  h  is  the  side  of  a  square  of  a  uniform  network,  so  that 
h±4  h(i  -  1,  2,  5,  4)  . 

In  the  case  of  boundary  conditions  of  the  general  type  it  is 
necessary  to  have  a  formula  to  approximate  the  normal  derivative. 

We  shall  derive  one  such  formula,  which  is  accurate  up  to  h3.  More 
accurate  formulas  may  be  obtained  by  using  relationships  (3l)»  (32), 
and  (^2)  for  two  spatial  coordinates. 

Let  us  consider  the  edge  of  a  plate  parallel  to  the  y-axis, 
and  let  us  set  up  a  relationship  for  the  boundary  condition  containing 
the  derivative  with  respect  to  the  direction  of  the  normal  perpendicu¬ 
lar  to  this  edge. 

For  this  purpose,  let  us  represent  u(x  +-h,  y,  t)  in  the  form 

,  -  ^  t  . 

u  ( v  jl.  y%  /)  =»  u  (.v,  r,  l)  huw  (.v,  r*  /)  -] (*»  JP»  0  +  Ry 

Expanding  the  functions  u(x,  y  —  h,  t)  and  u(x,  y  +  h,t)  according 
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to  Taylor's  formula,  we  find: 

*  (*.  y  +  4.  0  —  s»  (.v.  V, »)  +  u(.Y,  r  —  /«,  /)  -=  A*«',  (*,  y,  l)  -j-  Rt. 

Therefore,  If  we  discard  the  remainder  term  and  use  the  relationships 

thus  obtained,  we  may  write  the  equality 

S II  (.V  -f  h,  v,  0  -f-  II  (x,  V  —  A,  /)  —  411  (X,  r,  0  +  » (*,  y  —  A.  0™ 

“  2 fill',  (x,  Y,  /)  —  A*i«  /). 

Next,  let  us  write  out  the  same  equality  for  the  moment  t  +  l, 
add  it  to  the  preceding  one,  and  then  convert  the  relationship 
obtained,  using  for  this  purpose  Formula  (32)  written  for  two  spatial 
coordinates. 

As  a  result  we  obtain! 

.  (4  -r  2a)  u  (x,  y,  l  1)  m  —  2 h  \u't  (x,  y,  l )  'r  (x,  y,  t  ->  /)]  + 

-i-  “  (v»  y  -•  h,  t  -f  l)  -f.  211  (v  -M,  r,/4-  0  »  (v,  r  —  (  7l) 

(4  —  2^) 11  Cx*  y>  0  +  " (•''•  v  —  0  + «« (.v,  v  -r  A,  /)  -f-  2«  (x  -f-  A,  j,  /). 

If  Eq.  (42)  is  solved  with  the  aid  of  Formula  (70),  then  in  ( 71) 
it  is  necessary  to  take  X  =  4  and  to  replace  the  partial  derivatives 
with  respect  to  x  in  accordance  with  the  boundary  condition  on  the 
plate  edge  under  consideration. 

It  is  easy  to  derive  formulas  similar  to  (7l)  for  the  remaining 
three  edges  of  the  plate. 

C.  Let  us  consider  the  problem  of  the  propagation  of  heat  in  a 
body  bounded  by  a  surface  S,  when  the  initial  temperature  inside 
the  body  assumes  the  values  f(x,  y,  z) ,  while  the  boundary  (the 
surface  S  bounding  the  body)  is  maintained  at  the  temperature 
<p(x,  y,  z,  t)  at  all  t  >  0. 

It  is  required  to  calculate  the  temperature  of  the  body  at 
each  of  its  points  at  any  moment  of  time. 

The  solution  to  the  problem  reduces  to  integration  of  the 
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the  three-dimensional  equation 

(i.s4r+4r.+  *.\  (72) 

01  V  ox  0 y*  6?  J 

with  the  boundary  condition 

u  '■  ?  (xt  y»  t,  0  "Kw*  /  >  o 
and  the  initial  condition 

u  (.V,  y,  o)-f  (x,  V, 

Let  us  consider  a  cubic  network  in  a  space  with  a  rectangular 
coordinate  system  x,  £,  z  with  the  faces  of  the  cubes  of  the  network 
parallel  to  the  coordinate  planes. 

Let  us  now  take  a  cube  of  the  network  with  edges  2h  and  its 
center  in  the  node  marked  0  (main  node)  ;  its  upper  face  shall  be 
called  the  first  square;  the  cross  section  of  the  cube  cut  by  a  plane 
passing  through  the  point  0  parallel  to  the  upper  face,  the  second 
square;  and  finally  the  lower  face  shall  be  called  the  third  square. 

The  values  of  u  at  the  moment  _t  satisfy  [23]  in  the  network 
nodes  the  relationships! 


24«„»=2  H(, IT  "l,  |  —  6/lM  U„,t  -  0,5  —  R, 

i-i  * -7 

6  26 

>6h0.«=*‘s  H.-I+  y~'i  «„  I  — 12  AM  «o, I  —  A*  —  j?*. 


(73) 

(74) 


where  u0.t  is  the  value  of  u( x,  y,  z,  t)  in  the  main  node  O;  Z  u,  . 

i»l  x>z 

is  the  sum  of  the  values  of  u  at  the  mid-points  of  the  sides  of 
the  second  square  and  in  the  centers  of  the  second  and  third  squares; 


se 

2  1 
r=  19 


r,t 


denotes  the  sum  of  the  values  of  u  at  the  vertices  of  the 


is 


first  and  third  squares;  and,  finally,  Z  is  the  sum  of  the  values 

k-7 
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of  u  at  the  mid-points  of  the  sides  of  the  first  and  third  squares 
and  at  the  vertices  of  the  second  square.  The  errors  Incurred  from 
discarding  the  remainder  terms  in  (75)  and  (7*0  are  estimated  by 
the  inequalities! 


;pis 


21  i/«A* 


20 


where  Me  is  the  maximum  absolute  value  of  the  sixth-order  partial 
derivatives  of  u(x,  y,  z ,  t)  inside  S. 

Next,  we  use  (as  was  done  above  in  the  two-dimensional  case), 
together  with  ( 73)  and  ( 7*0  ,  the  Taylor  expansions  with  respect  to 


t  for 


** (*.  Jf.  1  l)  '?» (X,  y,  2.'), 


where  a  and  P  are  undetermined  factors,  and  we  choose  a  and  p  so  as 
to  obtain  relationships  which  are  accurate  to  he. 


We  thus  find  two  general  formulas.  The  first 

6 

.3«0,  i+li=  -  *"e,t+i  —  (24~ 

i  -l 

:8 

•f  /  *  8|l 


(75) 


where 

a-l2>.  — >.5,  y. 

for  which  we  shall  have 

jJKis;*:  -  >R»1  -* 

The  second 

3  M#,  *+« !<(,,  f-i-1  r  ^  Ui,t  —  y*  Wr, 
•  *  i—i  r-  19 

—(56— a • 
valid  for  any  a  and  P  such  that 


(76) 
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with  the  remainder  term  satisfying  the  Inequality 

./Si  ^ 

In  the  last  two  Inequalities  the  estimates 


l^ll  r~ 


9  31, A* 

2A3 


..-5 6 A/, A* 

t  •  — =  ”  «Ta  “  * 


should  be  fulfilled  for  Rj  and  Ra. 

When  X  =  6  there  result  as  a  special  case  from  Formulas 
(75)  and  (76)  the  formulas  examined  In  the  literature  [24]! 

6  iK 

*  r 

1 4 , 

s6 


Vk,  1-4-12 

p  s6 

£  u. '  • •• « + 16  r-  ]• 


The  coefficients  In  these  formulas  are  such  that  they  satisfy 

(§2)  all  the  requirements  for  the  convergence  of  the  computational 

process.  It  follows  from  this  same  Section  2  that  the  inequalities 

l-l  ~  - 1 

•  >0 

.  yrTU,*' 

!?l  :  *T - - • 

must  hold  true  for  the  errors  in  these  formulas  at  the  moment  T. 

Such  estimates  can  be  obtained  for  the  general  formulas,  if 
the  variations  in  X  are  limited  by  certain  inequalities  guaranteeing 
the  non-negativity  of  the  coefficients  in  the  formulas  under  considers 
tion. 

Thus  for  Formula  ( 75)  the  inequalities 

-»>o,  *  — 

should  be  fulfilled. 

Study  of  these  inequalities  shows  that  they  are  fulfilled 
simultaneously  for  each  X  of  the  interval 

!9-i- 1'jT- 

In  the  case  of  Formula  ( 76)  we  have 


D.  Here  also  the  problem  of  the  numerical  solution  of  the 
three-dimensional  heat-conduction  equation  may  be  reduced  to  a  system 
of  linear  algebraic  equations.  For  this,  we  take  Eq.  (73)»  discard 
the  remainder  term,  and  eliminate  A2U0jt  from  the  equality  obtained, 

using  for  this  purpose  the  relationship 

6 

h*  ]T  A  l\  t  -  6A>  A  l  ,  -p  £*.  (  77) 

1=1 


which  Is  obtained  from  Eq.  ( 12)  in  the  literature  [23],  If  we  replace 
u  by  A  u  In  ( 12) ,  multiply  the  equality  obtained  by  h2,  and  discard 
the  remainder  term. 

Then  we  can  write 


4» 


o 

s  ro,t- 4^ 


iX 

r„+.£  r 

'*7 


o 


»=i 


¥e  substitute  t  +  l  in  place  of  t^  In  this  formula,  and  the 
formula  obtained  to  the  original  formula,  replace  Au^  t(k  =  0,  1,...,6) 
by  ,  y^,  z^,  t)  In  the  sum,  and  use  an  equality  simi'  -  to 

(32)  written  for  three  spatial  coordinates. 

We  thus  obtain  the  relationship 


and  another 


(48+1 2k)  1 V  M4  -  (4  ~  2k)  Vi,  ,+l  -  2  Y]  r», ,+,  + 

•  -  I  fey 

6  ,8 

+ (4  +  =»>£  <-■+»  T  r,. ut-tW", 

*  4-  7 

'  6 

(56  4- 12  >•)  r0, ,+,-(»  -2X)  r<,  <+I + £  l+|  + 

o 

~  "r  2A^  Xj  l  i'  *  +  Xj  ,  r-*  ~  (>6  -  1 2 a)  l/*„  ,, 


r-19 


constructed  with  the  aid  of  (7*0  and  (77)  • 

The  number  of  these  formulas  may  easily  be  increased  by  using 
the  equalities  presented  In  the  literature  [23]. 
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For  example,  we  add  to  these  formulas  one  more  formula 

.  6  *  5#i 

< S24  —  8o X)  l\,  (4-1  —  5*  Vt>  >+»  -i*  (4  —  a '•)  ^  l\, *+4*-f 


r»i9 


(78) 


»n  u 

(4  +  apy**  l'r, ,  4-  ja  Ui, ,  —  (224  —  So a)U„  ,, 


rmjf 


1—  1 


hs 


which,  when  X  =  =  2,  assumes  the  form 

la2 


6  20  i> 

4^  U ,t I  +1  =  4  ^  ^11  (4-*  T  ^  Ur,  4  T“  4  ^  *  —  S  U%,  |.  (  79) 

.  1=1  r=  19  1=1, 

Formula  (78)  is  derived  in  the  same  way  as  the  preceding  two 
by  using  (7*0  and  the  equality 

2ft 

Ir'y  A  ,=  SA*Ar,»4  +  4^'^’ 

f  =  19 

which  we  find  using  Eq.  (17)  from  the  literature  [2J], 

Formula  (78)  leads  to  a  system  of  equations  which  may  be  s-lved 
by  the  method  of  successive  approximations.  Convergence  takes  place 
for  values  of  X  satisfying  the  inequality 

5 


We  arrive  at  this  inequality  by  following,  in  a  general  way,  the 
discussions  presented  above  for  the  two-dimensional  case. 

For  this  purpose,  let  us  consider  a  cube  with  faces  parallel 
to  the  coordinate  planes  and  located  at  a  distance  h  from  the  boundary 
under  the  assumption  that  Qj  lies  wholly  inside  the  given  cube  for 
which  Eq.  (72)  is  solved.  Network  nodes  lying  on  Qi  shall  be  called 
first-proximity  nodes.  Next,  we  take  a  cube  (lying  wholly  inside  Qt) 
again  with  faces  parallel  to  the  coordinate  planes  and  located  at 
a  distance  h  from  the  corresponding  faces  of  Qj,  and  the  nodes 
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lying  on  It  shall  be  called  second-proximity  nodes,  etc.  It  remains 
to  make  estimates  of  the  error  in  the  solution  along  the  cubes  $t» 
...»  on  which  the  nodes  lie,  in  exactly  the  same  way  as  we  estimated 
the  error  along  the  squares  7i,  7a,...  above. 

We  shall  now  solve  an  example  which  explains  how  the  method 
that  has  been  set  forth  is  applied  to  physical  problems.  Assume  that 
a  homogeneous  cube  of  cast  iron  with  side  L  *  1  m  is  being  cooled 
and  that  a  temperature  of  0°C  is  maintained  on  all  the  faces  of  the 
cube  throughout  the  entire  cooling  process.  It  Is  required  to  find 
the  temperature  distribution  inside  the  cube  at  the  moment  of 
time  T  =  1.6  hrs,  when  the  initial  temperature  is  distributed  (inside 
the  cube)  in  the  following  way! 

(80) 

u(.v  y,  o)  *  50  sin  ”.v  sm  sv  sin  ' 

In  order  to  solve  our  problem,  we  must  Integrate  Differential 
Equation  (72)  with  Initial  Condition  (80)  and  the  boundary  condition 

u  =  0  (81) 

for  t  0. 

In  determining  the  thermal  diffusivity  of  cast  iron  we  assume 
that  the  thermal  conductivity  is 


the  specific  heat 


the  specific  gravity 


a=54 


Veal 


hr  •  a  •  °C 


heal 


C,  a ltd 


r»72oo_* 

M 5 

where  kcal,  as  usual,  means  kilocalorie  ( large  calorie) . 
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Next,  let  us  proceed  to  calculate  the  thermal  diffusivity. 
According  to  the  data  of  the  preceding  paragraph,  we  find: 

X  j ti1 

m — -  — OtO<>2<  - . 

We  are  now  able  to  Integrate  Eq.  (72)  numerically.  For  this  purpose, 
let  us  take  the  length  of  an  edge  of  a  network  cube  h  =  0.2  and  use 
a  seven-term  equality  of  type  (79) ,  for  which  the  spacing  with 
respect  to  t  should  satisfy  the  condition 

J  =  — - — =0,32  (hour*)* 

2  a% 


The  Initial  values  of  the  temperature  in  the  network  nodes  are 


calculated  from  (80)  . 


i 


Fig.  2. 
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In  setting  up  the  equations  It  Is  advantageous  to  use  the  cross 
sections  of  the  cube  cut  by  the  planes 

— (i“  ti  ii  3«  4)> 

5 

perpendicular  to  Its  edges.  These  cross  sections  are  presented  in 
>  Fig.  2;  In  them  certain  nodes  [viz.  those  In  which  the  values  of 
u(x,  y,  z,  t)  coincide  due  to  the  symmetrical  temperature  distribu¬ 
tion]  are  designated  by  the  same  numbers. 

In  the  network  nodes  lying  on  the  boundaries  of  these  squares, 
as  well  as  on  the  faces  of  the  cube  z  =  0  and  z  =  1,  the  temperature 
at  any  moment  of  time  should  be  taken  equal  to  zero.  In  accordance 
with  (8l)  . 

Applying  Formula  ( 79)  to  our  example  ( to  the  nodes  designated 
in  Fig.  2) ,  we  obtain  a  system  of  linear  equations! 

4S  U 12U j,  <+!  —  SI  1,1  ~  12  U at  t 

44  cr.,i+i-4(r1,<44+3  r„  «+*)  ~aUi^  — 4^*.<+9  + 

40U',.|+|‘*4U\,|+|  rStr»,Mj)  +  9bT,,«-i-sOr3,*T- 

56  UA,  1+1=  12  I'l,  t+1  r  Ui, ,  +  )Ut,t  +  I  $  U}, I  -ir  t  , 

solving  which  for  t  =  0,  l,  21,  31,  and  4l,  we  find  a  solution  to 
Eq.  (72)  (with  coefficient  a2  =  O.O625)  such  that  when  t  *  0  It 
reverts  In  the  network  nodes  to  the  given  initial  state,  according 
to  ( 80)  . 

To  avoid  the  method  of  successive  approximations,  it  is  expedient 
to  first  solve  the  obtained  system  of  equations  for  the  unknowns 
Uk  *  1»  2,  3,  *0  ,  and  then  perform  the  computations  for  each 

layer  using  the  temperature  values  in  the  nodes  of  the  preceding 
layer.  We  obtain  the  following  four  formulas. 
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143^8  —  2121  f-rt.«  *r  35°4  t' j.i  -5*  849  Utu  4-  494'^«<( « 

I4j88  U»,4J**} l6S  l’i,i  —  372  ^V^"339<*  +  777  f,,i  , 

14388 U 283  r -j*  3  398  //,.« +  2067  Uy,  ~  2328  U^. 
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The  results  of  calculations  using  these  formulas  are  given  In 
Table  5. 

TABLE  3 

Temperature  values  In  nodes  of  cubic  network  calculated 
according  to  formula  ( 79) 

_  -  ■  T  *  j  : 

!  Appr<ud*at*  T*lu#»  '  Exact  vm lifts 


l 


I 

I 

IV  * 

. 

tj.i 

v, 

V„, 

U,.i 

tv« 

0 

10,15 

16,45 

2*1$* 

45,01 

IO.13 

1643 

26,58 

43»°i 

1 

5*55 

9,00 

14.^6 

25*55 

5.62 

9,09 

14,7° 

2>.79 

7l 

5*°4 

4  93 

7.97' 

12,90 

5.1* 

5.05 

8,13 

!3t*o 

.  */ 

1,07 

2,70 

4*57 

7*06 

1.72 

2,7& 

4.50 

7.28 

4/  ' 

0.91 

148 

-39 

3*87 

°»95 

1.54 

2Al 

4*°3 

5/. 

•  0,50 

0,81 

1.5* 

2,12 

0.53 

0,85 

0 

1.58 

2,23 

The  exact  values  of  u^^  In  the  node  (x^,  y^,  z^)  for  various 
moments  of  time  t  =  kl  are  calculated  using  the  formula 
Hi, *I*=  30 sills.- U sin sv,  sins y*»‘,,SM=*  (i=»o,  I,...',  j), 
since  the  function 

« (*.  v,  0  =  50  sin  ~  x sin  ~  y sin  n  -r-#’J,7i<  s* 

satisfies  the  given  differential  equation,  as  well  as  the  given 
Initial  and  boundary  conditions. 

The  method  Is  also  applicable  In  the  case  of  an  arbitrary 
closed  surface.  It  Is  only  necessary  to  set  up  special  equations 
for  the  boundary  nodes,  as  we  did  for  the  two-dimensional  case.  A 
large  number  of  such  equations  may  be  constructed.  For  example, 
the  simplest  of  them  can  be  obtained  with  the  aid  of  the  formula! 
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where  u„  u„  4.,...,  u-  4.  are  the  values  of  the  function  u(x,  y,  z,t) 

L|  v  * 

In  the  main  node  v  and  the  nodes  next  to  it  C,  A, . . . ,  E,  which  are 
separated  from  the  main  node  ( In  directions,  opposite  to  or  coincident 
with  the  directions  of  the  coordinate  axes)  by  a  distance  of  hj,  ha, 
...»  h6,  respectively;  and  Au„  .  Is  the  value  of  the  Laplace  opera- 
tlon  for  the  function  u  ¥  In  the  node  with  the  number  v  at  the 

V,  Ti 

moment  t. 

In  the  case  of  boundary  conditions  of  the  general  type  It  Is 
necessary  to  have  a  formula  of  increased  accuracy  for  the  approxima¬ 
tion  of  the  normal  derivative.  For  a  cube  they  can  be  obtained  by 
reasoning  in  the  same  way  as  in  item  B  in  the  derivation  of  the 
formula  for  a  square. 
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